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This  paper  reviews  attempts  at  signal  detection  in  Gaussian  noise  using  a  higher 
order  statistical  (Higher  Order  Spectra  (HOS)  or  polyspectra)  technique.  Examples 
comparing  power  spectral  and  blspectral  analysis  Include  the  following  topics: 
the  identification  of  signals  generated  by  a  system  of  coupled  nonlinear 
differential  equations,  radar  backscatter  processing  and  target  identification, 
and  a  statistical  treatment  of  the  detection  of  narrowband  harmonic  components 
resulting  in  a  Receiver  Operating  Characteristic  (ROC)  curve. 

The  critical  signal  and  noise  probability  density  function  (pdf)  assumptions  from 
polyspectra  theory  which  must  be  met  for  more  effective  noise  suppression  relative 
to  classical  second  order  power  spectral  methods  are:  (1)  discussed  in  relation  to 
detection  results  as  reported  in  the  literature  review;  and,  (2)  illustrated  via 
examples  using  both  direct  and  indirect  nonparametric  Discrete  Fourier  Transform 
(DFT)  bispectrums  employing  Fast  Fourier  Transforms  (FFT)  for  both  simulated  and 
real  data.  The  signal  set  consisted  of  pure  tones  and  a  hop  code  used  in  active 
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ABSTRACT 

This  paper  reviews  attempts  at  signal  detection  in  Gaussian  noise  using  a  hi^ier  order 
statistical  (Ifi^ier  Order  ^)ectra  (HOS)  or  ptdyspectra)  technique.  Examples  ounparing  power 
q)ectrBl  and  bispectral  analysis  include  the  following  topics:  the  identification  of  signals  generated 
by  a  sy^em  of  coupled  nonlinear  differential  equations,  radar  badcscatter  procesang  and  taiget 
identification,  and  a  statisdcai  treatment  of  the  detection  of  narrowband  hannonic  componems 
resulting  in  a  Receiver  Operating  Chancteristic  fflOQ  curve. 

The  critical  signal  and  noise  (xobability  density  function  (pdO  assumptions  ftor 
pcdyqKCtm  theory  which  must  be  met  for  more  effective  noise  suppression  relative  to  classical 
second  order  power  ^rectralmefiiods  are:  (1)  discussed  in  relation  to  detection  results  as  reported 
in  the  literature  review;  and,  (2)  illustrated  via  exanaples.  . using  both  direa  and  indirect 

nonperametric  Discrete  Fourier  Transform  (IXT)  biqpectrdriis  enq)loying  Fast  Fourier  Transforms 

» 

(FFT)  for  both  simulated  and  teal  data.  The  signal  set  Consisted  of  pure  tones  and  a  hop  code 

i  -  ^  ■ 

used  in  active  sonar.  [ _ 

The  results  are  in  agreement  with  the  examines  reviewed  whidi  concluded  that  no 
processing  gain  is  doived  fiom  bispectial  analysis.  It  is  shown  that  the  resolution  for  direa 
methods  uring  t»  formal  cumulant  coisttuction  (for  zero  mean  strictly  stationary  raidom 
processes,  cumulaius  up  to  third  order  are  etpul  to  moments  up  to  diitd  order)  was  far  stqretior 
to  foat  of  indirea  methods  containing  third  mder  cumulant  sequences.  Simuhtted  and  real  tests 
indicate  no  significant  improvement  in  biqrectral  over  power  qrectral  analysis  because  die  eOioed 
signals  exhibited  a  symmetric  pdf  (skewness  approximately  zero)  similar  to  Gaussian  noise.  In 
dieoiy  die  bispectrum  canna  differentiate  between  non-Gaussian  symmetric  signals  and  Gaussian 
or  non-Gaussian  symmetric  noise,  all  of  wtiidi  are  equally  suppressed.  Noise  ordy  suppression 


iv 


OCCUR  when  the  signals  echoed  or  radiated  are  zero-mean  non^laussian  with  non-zero  skewness 
(no  symmetry  evideitt),  and  die  noise  is  zero-mean  Gaussian  or  non-Gaussian  widi  zero  skewness 
(symmetry  evident). 

The  efficacy  of  signal  detection  in  underwater  acoustics  using  poly^iectral  methods  is 
examined  in  the  light  of  current  research  into  nonstadonary  HOS  with  Aiture  directions  for  study 
indicated. 
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CHAPTER  1 
INTRODUCTION 

Power  spectral  infonnation  is  that  wtuch  is  present  in  the  second-order  statistics  or 
autocorrelation  donuin  which  suffloes  fi>r  the  complete  statistical  description  of  a  Gaussian 
process  of  known  mean.  Phase  infonnation  is  suppressed  in  the  autoconelation  domain  since 
die  signal  in  qoestioo  is  repteneated  as  a  supeipositioo  of  stadsticaUy  uncorrelated  harmonic 
components,  ^tuations  may  exist,  however,  when  the  minimum-phase  assumptions  of  second- 
order  statistics  do  not  sufRce.  Le..  non-minimum  phase  infonnation  needs  to  be  extracted  as 
weD  as  deviations  fiom  Gaussianity  and  degrees  of  nonlineaiities.  Fm  these  cases,  hi^ier- 
order  qwctra  such  as  the  biqwctnim  and  trispectnun  are  required. 

Ifi^ier  order  spectre  or  polyspectre  are  the  (n-l)  dimensional  Fourier  transfoms  of  nth 
order  cumulant  sequenoes.^^  Ounulants  are  phue  smisitive  higher  order  statistics  which  are 
useful  in  die  analysis  and  desct^on  of  non-Gaussian  processes.  They  diqilay  die  degree  of 
higher  order  cnrelatiaa  present  in  a  random  series  as  well  as  providiiig  a  measure  of  the 
deviation  of  a  particular  distributkm  fiom  Gaussianiqf. 

In  genenl.  die  three  main  motivations  bddnd  die  use  of  ptdyspectre  in  signal 
processing  are:  (1)  Gausrian  noise  siqtpression  for  processes  of  unknown  qtectrum 
characteristics  in  detection,  parameter  estimadret,  and  dassification  prddems;  (2)  phase  and 
magnitude  response  reconstruction  of  signals;  and  (3)  detection  and  duracterization  of  time 
series  nonlinearities.  The  first  motivation  will  be  die  focus  of  this  paper  and  is  based  on  the 
property  that  for  Gaussian  processes  only,  all  Iti^ier-order  spectn  of  order  greater  dum  two  are 
idemically  zero.^^  This  leads  to  die  conjecture  dut  for  qjplications  vriiere  the  observed 


waveform  consists  of  a  non-Gaussian  signal  with  non-zero  dcewness  in  additive  Gaussian  noise 
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or  •  non-Gaitssian  signal  with  noc-zeio  skewness  in  non-Gaussian  wnse  with  zero  skewness, 
certain  advantages  may  exist  in  signal  d^ectitm  in  polyq)ectia  domains.^ 

Higher-Older  q)ectni  (polyq)ectia)  defined  in  tenns  of  the  higher-oider  moments  or 
cumulants  of  a  signal  (see  ^)pendix  B  for  definition  of  cumulants  and  joint  cumulants)  can  be 
used  to  examine  more  infonnation  contained  in  a  stochastic  non-Gaussum  or  d^eiministic 
signal  than  is  contained  in  its  autoconelation  (power  qectium).  From  statistical  dieory  we 
know  diat  all  joint  cumulants  higher  dum  the  second  order  vanidi  for  a  muldvariate  ncxmal 
process.^^  The  imtiMiHiate  impaa  of  this  is  diat  all  spectra  of  order  hi^ier  than  two  vanish  for 
a  Gaussian  process.  Thus,  higher-order  spectra  appear  useful  in  distinguishing  non-Gaussian 
processes  embedded  in  Gaussian  noise.  Moreover,  p(dyq)ectra  retain  phase  and  give  a 
measure  of  the  phase  cmrelation  between  components  whose  foequencies  sum  to  zero.  Such 
phase  relations  could  occur  due  to  nonlinearities  arising  from  modulation  effects  of  active 
sonar  signals  or  from  the  output  of  a  quadratic  system  diat  contains  power  contributions  at 
frequencies  which  are  sums  w  differences  of  pairs  of  input  frequendes,  a  phenomenon  known 
as  C^uadiatic  Phase  Coiqding  (QPQ.^^  QPC  could  result  from  the  interaction  of  diip  syrtems 
giving  rise  to  tonals  Oxnnonics)  above  broadband  noise  in  the  passive  sonar  case.  As  an 
example  of  this  consider  the  frequency  trijde  (fl4243sfl+f2).  Quadratic  oouiding  would  result 
in  a  power  contribution  at  f3  if  the  phase  at  f3  was  the  sum  of  the  phases  at  fl  and  fZ.  The 
biqiectium  identifies  quadratic  nrniUnearities  of  diis  type  triiile  the  power  qiectium  is 
completely  insensitive  to  such  phase  relatksis. 

There  are  basically  two  ways  that  can  be  used  to  estimate  lugher-order  qiectra  vriien  a 
finite  set  of  observation  measurements  is  available,  ctmventional  (Fourier  or  norqiarametric) 
methods  and  the  parametric  qiproach  based  on  Autoregressive  (AR),  Moving  Average  (MA). 
and  Autmegiessive  Moving  Average  (ARMA)  models.  Cmvemional  ^ifMoadies  can  be 
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imptememed  dtfaer  by  direct  or  indirect  mediods.  The  direct  method  refeis  to  direct  generation 
of  die  DPT  coefficients  required  for  die  ensemble  averse  of  the  bispectral  tii|de  (Moduct  from 
the  input  dtta.  It  does  not  require  uang  the  2*D  FFT.  The  indirect  method  refeis  to  the 
genenttion  of  cumulants  fiom  die  iiqMt  data  whidi  are  subsequently  double  Fourier 
trantfonned  to  give  die  bispectnim.  It  comes  from  the  general  definition  of  ptdyqiectia  being 
(n-l)  dimensional  Fourier  transfimns  of  nth  order  cumulam  fiinctimis.  In  the  case  of  bir^iectra. 
this  reduces  to  a  2-D  Fourier  transform  of  a  third  order  cumulant  function.  (See  sections  2.1.1 
and  2.12).  The  direct  mediod  as  a  ruUural  starting  point  for  polyspectra  analysis  is  die  main 
focus  of  diis  paper  widi  the  conventioiud  indirect  method  and  parametric  techniques  briefly 
examined  fiir  comparison  purposes. 

Few  qipUcations  of  polyspectra  have  arisen  due  to  the  voluminous  nature  of  the  data 
requirements  with  the  aocorrqianying  conqxitational  burden.  Also,  statistical  difficulties 
coupled  with  the  lack  of  a  physical  understandable  interpretation  have  compounded  the 
problem. 

For  our  purposes  we  will  only  be  concerned  widi  the  fact  that  if  a  time  series  has  a 
statistically  significant  bispectrum,  the  generating  process  is  non-Gaussian  with  a  suffident 
degree  of  skewness  as  to  be  prominent  in  a  Gaussian  or  non-Gaussian  symmetric  probabili^ 
dentiQr  function  Qxlf)  background  which  win  result  in  a  vanishing  biooherenoe.  (See  Appendix 
A).  Biooherenoe  analysis  results  when  the  bispectrum  shows  a  sensitivity  to  the  amiditudes  of 
die  involved  spectral  components.^  ft  has  been  shown  that,  for  even  inodierent  rignals,  die 
value  of  the  biqwctrum  cw  be  significant  if  averaging  is  perftHined  on  too  small  a  number  of 
records.  Since  the  bispectral  estimator  can  be  misleading,  normalized  bispectnims  or 
bioohereooes  are  the  preferred  choice  of  analysts  as  a  function  of  processing  ctqiabilities. 
Biooherenoe.  however,  is  very  sensitive  to  SNR  because  of  the  nramalization.  This  is 
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<4Merved  when,  for  perfectly  coherent  sipials,  the  bicoherence  is  not  one  due  to  the  input 
noise.  If  the  SNR  is  low,  the  value  of  the  bicoherence  estimated  with  a  relatively  small 
number  of  data  records,  e.g..  JiT  <  SO.  can  be  insignilicanL  To  recover  significant  value  the 
niunber  of  data  lecmds  must  be  increased  (up  to  several  hundred)  at  least  if  foe  signal  is 
stationary  over  a  sufficiently  long  time  interval.^  This  can  prove  prohibitive  which  is  the 
reascm  why  Inspectral  sensitivity  to  amplimdes  of  spectral  cmnponents  is  tolerated.  Moreover, 
udiat  is  an  undesirable  feature  for  some  is  advantageous  for  ofoers  as  this  same  amplitude 
sensitivity  is  foe  main  reason  that  uncoupled  harmonics  are  prominent  above  noise  in  foe 
bispectnim.  The  real  ptoUem  is  why  aren’t  such  signals  more  prominent  given  the 
"vanishing"  nature  of  bispectmms  and/or  bicoherences  to  Gaussian  processes.  The  literature 
would  seem  to  im{dy  that  a  processing  gain  is  inhetem  in  the  ^tfdication  of  polyspectral 
methods.  The  answer  is  that  theory  is  selttom  realized  in  practice  for  foe  reasons  cited  by 
Lagoutie  in  Reference  S  and  reported  here.  It  is  irttetesting  to  note  Lagoutte’s  conclusion, "  it 
can  be  said  foat  the  use  of  bispectnim  analysis  can  be  subjea  to  wrong  inteipretatioiL 
Biotfoerenoe  can  only  be  used  if  the  waves  are  not  embedded  in  noise."  This  author  has  not 
seen  too  many  physical  (xooesses  over  foe  course  of  his  career  foat  are  not  embedded  in 
significant  noise  of  some  type  other  that  those  idealized  results  which  invariably  turn  up  in 
journal  articles  but,  strangely  enough,  ate  never  applied  to  actual  {foysical  systems. 

*  Our  main  interest  in  the  use  of  polyspectra  is  the  realization  of  an  inqrrovement 
relative  to  second  order  methods  in  the  suppressicm  of  additive  Gaussian  noise  for  signal 
detection  problems  arising  in  underwater  acoustics.^  Non>Gaussian  signal  processes  occur  for 
active  sonar  when  the  reflecting  target  has  only  a  few  dominating  scatterers.  The  noise  in  sudi 
qrplications  is  fiequently  Gaussian,  so  foat  the  detectimi  proUem  is  foat  of  (toecting  a  non- 
Gaussian  signal  embedded  in  additive  Gaussian  noise. 
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Detection  {modems  related  to  non-Gaussian  noise^^  include  sonars  operating  under  ice 
udiere  noise  due  to  ice-cracking,  creaking  and  floe-smashing  contribute  a  compment  found  to 
have  substantial  non-Gaussian  behavior.  Additionally,  active  smats  operating  under  ice  near 
the  surface  may  encounter  a  mm-Gaussian  cmnpmient  due  to  qtecular  reflectitm  from  the 
irregular  under-ice  surface.  Another  environmental  situation  wdiich  may  produce  non-Gaussian 
noise  is  shallow-water  reveiberation.  Themetically,  non-Gaussian  noise  can  be  equally 
supfmssed  via  conventicmal  polyspectral  methods  if  the  noise  possesses  a  symmetric  pdf,  i.e., 
is  not  highly  skewed.^ 

This  pqrer  will  address  problems  related  to  active  sonar  as  defined  in  reference  22 
vriien  the  noiae  background  is  reverberation-limited  and  the  scatteiers  giving  rise  to  the 
reveiberatitm  are  assumed  to  be  statistically  independent  in  their  reflecting  pir^terties. 
^iplication  of  the  central  limit  theorem  gives  rise  to  a  Gaussian  process  for  reveiberation 
vdiile  the  return,  being  primarily  due  to  reflectimis  fit»n  a  few  random  scatteiers  will  be 
compositely  non-GaussiaiL 

In  passive  sonai^  die  detection  of  non-Gaussian  signals  in  Gaussian  noise  occurs 
when  the  background  noise  is  feequently  assumed  to  be  Gaussian  and  statitmaiy  with  the 
sources  such  as  ship-radiated  noise  being  non-Gaussiaa  The  pasrive  proUem  is  not 
considered  in  this  pqter. 

Ifi-Spec,  a  oollectitm  of  Matlab  M-files  designed  to  be  used  in  conjunctitm  with 
Matlab  developed  by  die  MathWorks,^^^^  is  a  software  package  provides  polyspectra 
analysis  ctqiabilities  for  various  i»ocessing  apfdicadcms.^’^  A  Mef  descr^ition  of  Hi-Spec 
routines  is  given  in  Appendix  A.  The  routines  ^lecifically  used  for  this  paper  included 
"bi^pec.d"  -  direct  method  for  Inspectnim  esdmaticxi.  "bispecj"  -  indirea  mediod  for 
biqwctiuai  estimation,  "gl.stat"  -  tfetection  statistics  for  Hiiddi’s  Gaussianity  and  linearity 
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tests,  "qpc_gen"  -  generator  for  quadratically  phase-cou|ded  harmonics  in  noise,  "tuum-gen"  • 
generator  for  harmonics  in  Gaussian  (colored)  noise  for  foe  hannmiic  retrieval  problem,  and 
"harm.est"  -  estimator  of  power  spectra  harmonics  usirrg  MUSIC,  Eigenvector.  Pisarenko,  ML 
(Cqxm)  and  AR  methods  based  on  the  diagonal  slice  of  the  fourth  order  cumulant  with 
ctmventional  periodogram  also  give  for  comparison  purposes. 

The  routine  "qpc_gen"  generated  foe  data  for  Figures  C.1>1  and  C.1-2  in  Appendix  C 
as  described  in  section  C.1  DIRECT/INDIRECT  COMPARISON  FOR  SIMULATED 
DATA. 

The  routine  "biqxc.d"  was  used  to  calculate  the  direa  bispectrum  for:  (1)  Figures 
CO-l,  C.0-2,  and  C.0-3  in  Appendix  C  as  described  in  section  Ci)  POWER 
SPECTRUM/BISPECTRUM  COMPARISON  FOR  SIMULATED  DATA.  (2)  Figures  C.l- 
1  and  Cl-2  in  Appendix  C  as  described  in  section  Cl  DIRECT/INDIRECT  BISPECTRUM 
COMPARISON  FOR  SIMULATED  DATA,  and  (3)  Figures  C2-1  and  C2-2  in  Appendix  C 
as  described  in  section  C  J  POWER  SPECTRUM/BISPECTRUM  COMPARISON  FOR 
REAL  DATA. 

The  routine  "bi^iecj"  was  used  to  calculate  the  indirect  Inspectrum  for  Figures  C.1-1 
and  Cl-2  in  ^ipendix  C  as  described  in  section  C.1  DIRECT/INDIRECT  COMPARISON 
FOR  SIMULATED  DATA. 

The  routine  "gl.stat”  was  used  as  an  added,  albeit  redundant  chedt  due  to  foe  observed 
obvious  non-zero  nature  of  the  edioed  signal  skewness,  to  verify  the  zero  skewness  hypothesis 
for  the  real  data  as  described  in  section  CJ2  POWER  SPECTRUM/BISPECTRUM 
COMPARISON  FOR  REAL  DATA.  The  use  of  the  "gl_stat"  routine  in  fois  context  was 
more  for  security  than  an  absdute  requiremem  and,  as  such,  not  formally  addressed  in  this 
paper.  It  will  suffice  to  say  that  "gl_stat"  verified  the  zero  skewness  hypothesis  in  all 
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instances  which  gave  die  author  a  warm  feeling  about  the  "correctness"  of  the  Hi-Spec 
software  package.  It  is  worth  noting  that  there  are  many  statistical  tests  of  a  less  craiplicated 
nature  then  "gl_stat"  vdiidi  serve  the  same  purpose,  i.e..  the  comparison  of  continuous  or 
Unned  data  where  one  data  set  is  compared  to  a  known  distribution  or  two  equally  unknown 
data  sets  ate  to  be  compared.  See  Qu-Square  goodness-of>fit  and  Kohnogorov-Sminurv  tests 
as  described  in  reference  52. 

The  routine  "harm_gen"  was  used  to  generate  fee  harmonics  for  fee  higher  order 
statistical  resolution  improvement  routine  "harm.est"  as  repotted  in  section  22.1  Harmonic 
Retrieval.  The  results  were  taiefly  addressed  for  iUustrative  purposes  to  make  fee  reader 
aware  that  higher  order  statistical  versions  of  these  routines  suffer  fiom  fee  same  constraints 
that  feeir  lower  order  counterparts  do  >  namely  a  high  SNR  requirement.  This  coupled  wife  a 
lack  of  signal  skewness  caused  no  discemaUe  difference  wife  second  order  periodogram 
methods. 

A  condensed  version  of  feis  pqrer  was  puUished  in  The  U.  S.  Navy  Journal  of 
Underwater  Acoustics,  Volume  43,  No.  1,  January  1993,  pp.  201-220,  under  the  title,  "Signal 
Detection  Using  fee  Biqxctrum." 
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CHAPTER! 

A  REVIEW  OF  DETECTION  RESULTS  AS  REPORTED  IN  LITERATURE 


This  chapter  reviews  noiq>anmetric  and  parametric  attempts  at  signal  detection  in 
Gaussian  noise  using  the  tM^rectium.  Examides  horn  the  literature  whidi  address  comparison 
of  power  qrectial  to  bispectial  mediods  include  the  following  q)plications:  (1)  die 
identification  of  signals  generated  by  a  system  of  coiqded  nonlinear  differential  etpiations,  (2) 
radar  hackscatter  processitig  and  target  identification,  and  (3)  a  statistical  treatment  of  die 
deiecdon  of  narrowband  harmonic  con^xments  resulting  in  a  Receiver  Operating  Characteristic 
(ROQ  curve.  The  nonparame^c  direct  and  indirect  methods  are  formally  defined  along  widi 
cofTC^XMidirig  test  results  sridch  are  in  agreemett  with  die  examides  cited  in  the  literature 
review.  Pumnetric  mediods  are  described  including  the  reasons  for  ne^ecdng  same  in  this 
p»ptT.  Harmonic  retrieval  is  briefiy  discussed  in  the  context  of  conjectured  resoludon 
enhancements  doe  to  qiidying  higher  cK&a  statistics  to  standard  algorithms.  The  resoludon 
section  is  included  to  show  the  extension  of  the  detection  pnddem.  ladt  of  signal  ricewness 
cooqiared  to  noise,  to  die  resoludon  protdem. 

2UU  Idariification  of  Signals  Generated  by  a  Coupled  System  of 

Ndnliiiear  Dlfllerential  Equations 

The  detectability  of  signals  was  examined  dirou^  eiqieriments  in  which  simulated 
signals  were  subjected  to  power-qiectral  and  direa  bispectial  analysis.^*  ^lecifically.  rignals 
lidi  in  harmonic  content  were  generated  by  solving  a  system  of  coiqded  ntmlinear  diffeiendal 
equations  and  were  subaequendy  contamiiuted  by  Gausrian  white  noise  before  being  subjected 
to  tpcanl  and  Uspectral  analysis.  The  eHects  of  start  iq>  transients  due  to  initial  conditions 
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were  removed  by  diiicarding  the  fiist  2JKO  poims  as  the  scdution  was  stepped  for  typically 
10,000  stqM.  Classical  windowing  and  segment  averaging  was  apf^ed  to  the  noise  cases  for 
two  different  noise  levels  described  as  low  and  hi^  In  both  instmces  ensemUe  averaging 
did  not  produce  sufficient  reduction  of  background  noise  to  make  the  bispectnim  a 
significantly  mtHe  efficient  signal  d^ector  than  the  power  spectnun  It  was  concluded  that  die 
biqiectnim  slices  taken  did  not  show  it  to  be  a  very  valuaUe  indicator  of  die  existence  of 
multiple  hannonic  signals  in  the  presence  of  Gaussian  noise.  No  signals  were  observed  in  the 
presence  of  noise  using  Mspectra  when  they  were  not  also  visible  in  the  power  spectra.  It  was 
further  shown  that  no  inoease  in  processing  gain  is  derived  from  bispectral  arulysis  even 
though  it  makes  use  of  inter-fiequency  phase  relation  data.  It  was  found  that  reliable  detection 
tequired  m<He  signal  energy  in  die  bispectral  case  than  for  qiectral  detection  unless  the  signal 
skewness  was  large.  The  amount  of  avoaging  applied  to  reduce  random  sami^ing  variations 
is  pngioitional  to  the  time  bandwiddi  product  The  disadvantage  of  bispectral  arulysis  varied 
direcdy  with  die  time  bandwidth  product  as  results  indicated  diat  die  bispectral  detector  was 
inferior  to  both  energy  and  correlation  (matdied  filler)  detect(»s.  Bispectra  did  not  appear  to 
offer  any  advantage  in  detectability  at  low  signal-to-noise  ratios  urdess  die  signal  skewness 
was  large. 

As  a  caveat  it  was  pttinted  out^*  that  the  signal  biqiecttum  contained  more  structure 
than  die  power  qiectrum,  indicating  that  bispectral  arulysis  of  signals  may  provide  useful 
additional  dassificatitm  beyond  that  obtained  frmn  enmgy  mediods.  It  was  further  conjectured 
that  biqwctral  analyris  nuy  be  useful  in  active  acoustic  detectors  in  the  presence  of 
reverberation  assutning  that  die  acoustic  return  reflects  a  localization  conqured  with  the 
reflectors  producing  leverbeiation. 
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Radar  Baducatto’  Processing  and  Target  Identification 
The  use  of  bispectral  processing  methods  for  radar  backscaner  ptocessii^  and  taiget 
identification  have  also  been  examined.  In  this  situation  bispectral  processing  methods  were 
adapted  to  the  radar  signature  analysis  proUem  which  resulted  in  the  birange  profile,  a  2-D 
display  of  target  scattering  mechanisms  in  the  range  domain  widi  the  final  goal  being  die 
detection  and  dasdfication  of  radar  signals. 

The  type  of  signal  processing  addressed  in  this  study  is  based  on  a  spedfic  target 
scattering  model  idiere  the  scattered  signal  is  a  comtxnation  of  specular  scattering  terms  from 
localized  areas  on  the  target  and  multiple  interaction  terms.  Ifigher  order  qiectral  processing 
of  radar  signatures,  similar  to  dassical  spectral  analysis  of  radar  data,  produces  signatures  in 
the  time  domain  that  ate  related  to  the  gemnetrical  shape  of  the  target  This  serves  as  an  aid 
in  target  discrimination.  A  description  of  this  process  follows. 

Complex  natural  resonances  of  radar  targets  are  aqiect  independem  and  are  used  as  an 
a^iect  insensitive  method  far  discriminating  between  targets.  No  information  about  the 
scattering  medianisms  of  a  radar  target  is  directly  evident  from  the  backsca tiered  frequency 
reqxmse  data.  However,  this  informatitm  pertaining  to  target  sluqie,  size  and  orientation  is 
present  in  fiequency  and  can  be  extracted  1^  exsnining  the  target  scattering  mechanisms  that 
are  disidayed  in  the  Target  Inqwlse  Re^xmse  (TIR).  The  TIR  comes  into  play  since  if  a  i^ane 
wave  is  transmitted  to  illuminate  a  target,  then  the  backsca  tie  red  field  is  a  fiinction  of  the 
transfer  function  as  seen  the  radar,  and  the  range  to  die  target  Accordingly,  it  is  desinMe 
to  change  variables  fimn  frequency  to  time  since  target  scattering  features  can  be  recovered 
fiom  the  impulse  respemse  giving  the  analyst  a  more  intuitive  understanding  of  the  processes 
invrdved  by  working  with  time  domain  signatures  instead  of  frequency  domain  reqxmses.  For 
the  radar  proUem  the  end  remit  beemnes  a  time  trqde  correlation  mqiped  to  a  range  idane  via 
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die  biqwctrum  where  the  bispectnim  is  expressed  ss  an  ensemble  average  of  impulse  responses 
B(rl  j2)  «  <di(rl)Ji(i2)Ji(rl-i-i2>  where  h(r)  is  die  impulse  response  as  a  funcdon  of  range. 
Defining  the  Urange  as  such  implies  that  a  non-zero  biqiectnl  response  is  die  result  of 
interaction  between  the  responses  at  ranges  rl  and  i2  udiich  appear  as  a  tespoost  at  the  range 
rUt2.  Fbr  the  (nqiectnim  of  real  data,  this  usually  takes  the  fonn  of  a  cleariy  evident  non¬ 
zero  magnitude  bispectnim  at  die  ordered  pair  in  range  (rl  j2)  in  the  principal  triangular  region 
of  die  tn^iectnun  (see  qipendix  B),  and  various  assorted  other  ordered  pairs  in  the 
coireqionding  eleven  biqiecttum  symmetry  regions.  (See  Figure  B.2-lb.)  These  symmetry 
regions  containing  significant  non-zero  magnitude  biiqiectrums  occur  at  the  ordered  pairs: 
(r2jl)  from  symmetry  regim  2.  (-i2,rl-»-t2)  fiom  conjugate  symmetry  region  3.  (-rl/l-ft2) 
from  oonjngate  symmetry  region  4.  (-rl-r24'l)  from  symmetry  r^on  S,  (-rl-t2,r2)  from 
symmetry  r^on  6,  (rl.-t2)  from  conjugate  symmetry  region  7,  (-i2,-rl)  from  conjugate 
symmetry  region  8.  (r2,-ri-i2)  finom  symm^  region  9,  (rl,-rl-r2)  from  symmetry  region  10. 
(rl-fr2,-rl)  fiom  conjugate  symmetry  region  11.  and  finally.  (rl-i-i2,-r2)  fiom  conjugate 
symmetry  region  12.  (See  Appendix  C  section  Cl  figures  for  an  illustration  of  the  twelve 
symmetry  regimis  of  the  bispectrum  for  real  data.)  Furthermore,  die  bispectrum  at  (rl/2)  is 
mm-zero  only  if  the  re^xmses  at  rl  and  r2  and  rl-t-t2  are  correlated.  For  die  biqiecttum  of 
comfdex  data,  the  non-zero  magnitude  of  the  biqiectral  reqxmse  coirespmiding  to  die 
interactions  between  the  reqxmses  at  ranges  rl  and  i2  sdll  qipears  at  the  ordered  pair  in  range 
(rl  j2)  located  in  the  principal  triangular  r^on  of  die  biqiectrum  as  before  for  the  real  data 
case.  However,  because  we  are  now  dealing  with  bispectrums  of  complex  as  t^iposed  to  real 
data,  there  are  no  longer  twelve  symmetry  regimis  but  only  one  defined  by  die  idacement  of 
die  comidex  conjugate  for  die  third  order  jmnt  cumulanL^*^^  An  exaoqde  of  die  idacement 


used  in  Hi-spec  is  given  in  ^ipendix  B  diiecdy  below  Rgure  B.2-3.  This  corresponds  to  the 
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dm  symmetry  relation  dedned  ft>r  die  bispectium  of  complex-valued  signals  immediaiely 
following  the  exanqde.  Foimally,  this  range  ^mmetiy  fx  dre  complex  case  would  only  occur 
at  (i2jl)  in  the  area  designated  region  2  in  Rgure  B.2-lb  in  ^ipendix  B.  The  noise 
siqipression  ci^Mbilides  of  die  biqiectram  arise  ftom  the  firet  that  unomrelated  or  low 
cmrelated  Gaussian  noise  of  unknown  qiectra  diarsctenstics  (coloied  or  white)  does  not  have 
a  signidcant  bispectrum.  Gonaequendy.  the  omtivation  for  using  bispectrums  for  detection 
purposes  becomes  dear  as.  in  theory  for  an  indnite  numbo*  of  samples.  Gaussian  processes 
varddL  hi  practice,  however,  only  a  finite  number  of  data  samples  is  available  so  additive 
Gaussian  noise  is  not  totally  suppressed  in  the  binnge  profile.  Also,  since  additive  noise  may 
be  only  asymptotically  Gaussian  or  non-Gaussian  dutter  having  a  non-aeio  dtird  order  momem 
(dtewness).  total  noiae  suppression  wouldn’t  occur  for  even  an  indnite  number  of  samples. 

The  obvious  conjecniie  is  that  the  bi^KCtnan  miilit  prove  a  viable  ahemative  to  power 
yectrum  tedutiques  since  it  retains  phase  infonnation  thereby  identifying  any  type  of  phase 
coupling  udiidt  might  result  dom  the  interactions  of  target  scattering  mechanisms.  This 
information  could  allow  ftv  better  target  detecdoiL  The  time/range  mapping  is  discussed  next 
In  most  dgnal  processing  problems,  die  data  sequence  consists  of  samides  taken  ftom 
a  time  dqiendent  waveform.  Then  the  tdqiectium.  as  dedned  in  ^ipendix  B.  is  a  functitm  of 
two  variables  in  die  fiequency  domain.  The  radar  scattering  data  sequence,  however,  is  not  a 
time  aeries  but  recorded  in  the  frequency  domain.  The  tr^de  correlation  of  diis  radar 
frerpiency  domain  data  is  a  2-D  prodle  in  die  frequency  donutin  with  the  qtfdied  biqiectium 
becoming  a  prodle  of  target  scattering  dgnatures  in  the  2-D  time  domain  whidi  could  be 
termed  a  "bi-time”  prodle  of  radar  targets.  This  prodle  can  be  subsequendy  eiquessed  ma2- 
D  prodle  in  range  using  die  standard  range  equation  r  B  ct/2  udiete  r  is  die  range  dom  die 
radar  to  the  target,  t  is  die  time  needed  for  the  signal  to  propagate  to  the  target  and  bade,  and  c 
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if  die  speed  of  light  This  is  is  lefened  to  u  the  binnge  profile  of  the  target  Formally, 
it  is  a  bispectral  display  oi  target  signatures  in  the  2-D  ran^  domain  rl  and  i2  where  die  term 
"range"  denotes  the  propagation  distance  to  the  target  This  process  is  direcdy  related  to  the 
classical  problem  of  determining  the  system  inqailae  response  by  taking  the  inverse  Fourier 
transfonn  of  the  system  transfer  funcdoo  (arudogous  to  die  radar  system  transfer  funcdon 
containing  the  complex  natural  resonances  of  die  targets).  This  resulting  time  triide  correlation 
is  subaequendy  mapped  to  range  to  create  die  birange  profile. 

The  focus  of  the  report  was  on  the  processing  of  radar  data  treated  as  a  time  series 
using  higher  Older  rpectral  arudysis  tedmi^ies,  in  particular,  the  biqrectrum.  The  reason  for 
the  renewed  interest  in  the  bispectrum  as  a  general  signal  processing  technique  is  reimed  to 
current  increases  in  processing  capabilities  as.  previously,  excessive  computadonal 
requircmeots  gave  the  promise  of  limited  rewards.  Not  only  are  die  madiemadcs  of  the 
bi^Kctium  for  more  comidicated  dun  for  mher  spectral  arudysis  techniques  but  the  physical 
insight  into  die  technique  is  not  fiiUy  understood  as  the  main  problem  remains  the  lade  of  an 
intuhive  interpretation  and  understanding  of  die  biqrec&al  processing  of  a  time  series. 

Both  Diiea  and  Indirect  mediods  of  biqiectium  computation  were  outlined  widi  die 
requhement  fiv  a  large  number  of  data  samides  due  to  the  hi^  variance  incurred  by  both 
estimates,  It  was  noted  diat  die  esdmates  could  be  made  smoother  if  mme  data  segments  were 
used  at  die  expense  of  decreased  biiange  lesdudon  diat  results  from  die  introduction  of 
nonstadonary  problems.  The  main  advantage  of  these  methods  is  the  ease  of  implemeiaadon 
using  HT  algorithms.  A  vaiiam  of  these  methods  was  used  to  generate  die  birange  profiles. 

A  parametric  technkpie  using  autoregressive  modeling  was  introduced  as  an  altemadve  to 
classical  Fourier  transfonn  tedmiques.  The  reasons  for  this  stemmed  from  the  resdudon  of 
the  estimated  profiles  being  limited  by  bandwidth  and  die  effects  of  die  window  function  used. 
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(The  window  had  a  wide  mainlobe  which  fimher  limited  reacdmkm.)  It  was  pointed  out  dut 
applying  a  rectangular  window  to  the  triple  correUttion  for  the  biqrectnim  2-D  skklobes 
to  appear  along  each  mge  axis.  It  was  obsenred  that  with  finiie  sets  of  radar  scattering  data, 
rescriution  pioblenis  widi  Fourier  tedmi^ies  are  inevitable.  To  avoid  this  problem  the 
parametric  technique  was  imposed  which  assumed  that  the  scattering  as  a  function  of 
frequency  satisfies  a  model  whose  parameters  are  estimated  from  the  measured  data.  The 
biruige  teairiution  was  no  longer  limited  by  the  measurement  bandwiddi  as  was  die  case  with 
Fourier  baaed  processing  of  finite  length  data  reomls.  The  parameters  in  (piestion  were 
autoregressive  (AR)  and  could  be  estimated  fiom  the  second  or  ddrd  order  cumulsnts  of  the 
measured  backacatter  data  sequence.  The  birange  was  then  computed  as  die  triple  corrdatkm 
of  the  inverse  Z  transfixm  of  the  measured  badcscatter  r^xnse  which  became  the  Target 
Impulse  Response  (HR).  Hence,  the  birange  was  parameterized  by  die  AR  coefficients.  The 
fimitatfoas  of  parametric  modeling  arose  fiom  die  faa  that:  (1)  scattering  u  fiequencies 
outside  die  measurements  can  not  be  easily  piedicied  from  the  measured  dam  due  to 
dependence  on  msiqr  fiKtors  such  as  dispn^  and  scattering  re^on.  aid  (2)  modd 
parameters  are  often  estimated  using  nonlinear  algorithms  and,  as  sudi,  are  sensidve  to  noise 
in  the  data.  It  was  noted  diat  improvements  in  Fourier  methods  could  be  addeved  if  die 
measured  hackscatter  had  a  high  signal-to-noise  ratio  coupled  widi  a  slowly  varying  fiequency 
leqxnse  fiw  the  target 

An  AR  parametric  model  where  the  parameters  were  estimated  using  third  order 
cumulants  was  devised  and  tested  using  experimental  radar  data.  The  AR  birange  profiles 
were  oooqiated  to  die  Fourier  transfimh  based  birange  estimates.  The  results  obtained  on 
modding  scaled  modd  aircraft  using  die  AR-based  algorithm  were  not  encouraging  in  the 
sense  that  no  dosdy  tptced  peaks  in  the  birange  were  resdved  and  little  compatibility 


15 


between  Foorier^tued  bumge  profiles  and  AR-teaed  binoge  imfUes  was  otoerved.  This 
incompatibUity  was  attributed  to  problems  associated  witb  modeling  complex  (bispectial,  ix., 
binnge)  as  opposed  to  spectial  data.  For  a  detailed  discussion  on  the  problems  associated 
with  the  AR  modding  oi  oom{dex-valued  signals  see  sectioo  4.3  of  Jouny.^ 

The  motivation  for  biqtectral  processing  of  tadar  signals  for  this  smdy^  was  the 
inherent  advantage  over  spectral  processing  techniques  for  the  suppressing  of  additive  Gaussum 
noise  or  any  additive  undesirable  random  signal  wifo  a  symmetric  pdf.  This  engendered  the 
interest  in  using  biiange  profiles  for  target  recognition  puipoees,  in  paiticular,  if  the  backscatter 
wavefimns  are  oomqtted  with  zero  mean  additive  Gaussian  noise.  The  process  of  target 
recognition  was  appioadied  in  two  ways:  magMuametric  FFT-baaed  and  parametric  AR 
modeled.  The  parametric  qtproach  required  the  knoudedge  of  die  underlying  distribution  of 
die  biiange  profile  and  the  apiiori  piobabilhy  of  occurrence  of  eadi  target  Such  an  approach 
was  deemed  useful  for  estimating  die  optimal  classification  peifi»inance  that  could  be  achieved 
using  the  biiange  as  a  feanne.  However,  due  to  the  lack  of  knowledge  of  the  apriori 
piobabilides  of  posriUe  targets,  and  due  to  the  computadonal  burden  imposed  by  parametric 
classification,  nonparametric  dassifiers  are  often  more  appropriate. 

In  summary,  it  was  found  that  peifonnance  gains  using  die  biiange  profile  were  slight 
with  the  addidonal  computadonal  needs  being  m  undesirable  tradeoff.^ 

2A3  Statistical  Treatment  of  the  Detection  of  Narrowband  Harmonic 

Components  Resulting  in  a  Receiver  Operating  Characteristic  (ROC)  Curve 
A  comparison  of  die  detecdon  performance  of  the  power  spectrum  and  bispectium  for 
delecting  narrowband  harmonic  components  is  given  in  IKfilsoiL^  The  onphasis  here  is  on  a 
more  sadstical  treatment  of  die  detection  process.  Consislent  estimates  of  the  urmoimalized 
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bispectrum  are  omstnicted  as  a  function  of  the  number  of  IXTs  over  wbich  the  luspectnun  is 
cdierently  averaged.  Nonnalized  bi^ctiums,  ix.  Incoherences  [see  (A.2-1)  and  B.2*l0)l.  are 
also  considered.  The  variance  of  the  unnoimaiized  lMq)ectial  estimator  bectnnes  a  function  of 
the  multiple  of  the  power  plectrums  at  the  respective  Inqsectral  frequencies  f(j),  f(k),  and 
f(j4k)  scaled  by  the  DFT  size  divided  by  the  number  of  DFTs  averaged.  The  time  series  is 
assumed  to  be  uncorrelated  in  time.  The  coherently  averaged  estimator  of  the  iMqrectrum  is 
shown  to  be  an  element  of  a  multivariate  comidex  Wishart  distribution  with  dimension  3, 
degrees  of  freedom  corresponding  to  the  number  of  averaged  bispectxums.  and  ctnnplex 
Gaussian  distributed  in  the  asymptodc  sense.  The  bicoherence  is  dttermined  to  produce  a 
quantity  whose  asymptotic  statistics  are  more  easily  calculated.  It  is  pointed  out  that  since  the 
a^mptotic  (Ustributitm  of  the  bi^ctral  estimator  is  cranptot  Gaussian,  the  distribution  of  the 
bictdietenoe  is  given  by  an  asymptotically  noncentral  clu-square  distribution  with  two  degrees 
of  fieedom  and  a  nrmcentrality  parameter  wfaidi  is  a  function  of  the  dcewness  of  tire  dau  in 
the  time  series  in  question.  e.g.  if  the  time  series  is  Gaussian,  die  skewness  function  will  be 
zero  for  all  biqiecttal  frequency  pairs  and  the  distribution  of  die  tdcoherence  (normalized 
bispectrum)  is  central  chi-square  with  two  degrees  of  fineedtnn  (the  noncentiality  paramet^  will 
also  be  zero).  For  the  detection  of  signals  in  the  presence  of  additive  Gaussian  noise,  a 
threshold  is  qiplied  to  the  tdadwrenoe  estimate  to  detect  a  non-zero  Uqiectrum  at  a  specified 
false  alarm  rate  ootreqxmding  to  the  probability  of  ctmimitting  a  Type  I  error  (accepting  die 
alternative  hypothesis  when  the  null  is  true)  from  die  central  dii-square  distribution.  The 
derivation  of  the  noncentrality  parameter  as  a  function  of  the  skewness  of  the  signal  and  the 
signal-to-noise  power  ratio  for  a  noncentral  chi-square  distributed  nonnalized  tn^pectrum 
estimate  of  signal-[dus-noise  is  given  in  Hinich.^^  The  probability  of  detectkm  is  then 


cmnputed  for  the  normalized  bispectrum  (Idctdierence)  from  die  noncoitral  dii-square 
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distribution  as  a  fimction  of  the  noocentrality  and  the  threshtdd  paiameteis  for  the  noncential 
chi-square  distribioion.  A  similar  threshold  is  derived  as  a  function  of  die  signal-to-noise 
power  ratios  at  the  f(j).  f(k).  and  fCH-k)  biqxctral  frequencies  for  die  magnitude  squared 
uimoimalized  bispectrum  to  calculate  the  probability  of  detection  from  the  noncential  chi- 
square  distribution  by  evaluating  the  probability  that  a  noncentral  du-square  random  vaiikUe 
with  a  given  non-centrality  parameter  will  exceed  the  normalized  bispectium  threshold  which 
is,  in  turn,  normalized  1^  the  aforementioned  bispectral  fiequency  signal-to-power  ratios. 
Stanley  put,  the  coqcon  is  how  large  die  normalized  biqiectnim  statistic  has  to  be  in  order  to 
confidendy  reject  the  "Gaussian  noise  only"  hypothesis  in  favor  of  declaring  the  presence  of  a 
non-Gaussian  signal  It  is  shown^  that  to  operate  at  a  false  alarm  rate  of  .001,  one  would 
reject  the  Iqrpodiesis  that  only  noise  is  present  (the  normalized  biqiectnim  statistic  is  central 
chi-aquare  distributed  radier  than  nonoentral  du-stpiare  distributed)  fm  values  die 
normalized  biqiectrum  statistic  that  are  13.8  ot  larger.  The  observation  is  made  diat  since  the 
mean  value  of  die  normalized  bi^Kctrum  statistic  is  equal  to  the  noncentrality  parmneter,  it 
follows  that  the  noncentrality  parameter  must  be  13.8  or  larger  to  addeve  detection  at  this 
false  alarm  rate.  The  critical  issues  addressed^  hom  a  statistical  perqiective  pertain  to  this 
noncentrality  parameter,  in  particular,  the  several  factors  which  contribute  to  its  value- 
skewness  (diaracteiistic  of  the  signal),  signal-toHioise  ratio  (characteristic  of  the  signal  and 
noise  power  levels),  and  finally,  the  i»ocessing  characteristics,  ix.,  DFT  size,  number  of  DFTs 
used  in  biqiectrum  averaging. 

The  oondusions  reached^  can  be  summarized  as  follows.  If  a  signal  has  a  mm-zero 
skewness  function,  a  tradeoff  exists  between  signal-to-noise  ratio  and  processing  parameters 
wUdi  will  determine  if  the  non-zero  skewness  is  nifficiendy  large  to  result  in  a  noncentrality 
parameter  which  wOl  aUow  its  detectkm  at  a  given  ftise  alarm  rate.  For  small  values  of  the 
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skewness  fiinctkn,  tbe  pfocessing  psrameten  snd  signal-toHn^  ratios  have  to  be  such  that 
their  pnxhict  remains  laigie  enough  to  prodi^  a  sufficiently  huge  ooncentrality  parameter  for 
detection.  [The  noncentrality  parameter  has  a  linear  dqtendence  on  the  mimber  of  averages 
and  an  approximately  cubic  dqtendence  on  signai-to-raMse  ratio  (for  low  Shot).]  The 
imfdication  here  is  that  if  the  signal-to-noise  ratio  decreases  by  a  factor  of  2  (3  dB),  then  it  is 
necessaiy  to  mcrease  the  number  of  averages  by  a  factor  of  8  in  order  to  retain  the  same  level 
of  detecttritility  dictating  a  certain  degree  of  flexibility  for  real  time  algoiidunic  realizations 
which  may  not  be  feanUe.  By  far.  it  is  clear  that  the  underlying  concern  remains  signal 
skewness  as  IKfllson^  demonstrates  the  essential  relationship  between  signal  characteristics, 
noise  characteristics,  and  processing  parmneteis  wtudi  define  the  detection  performance  of  the 
bispectrum.  It  is  shown  that  given  a  skewness  function,  the  processing  parameters  necessary 
to  addeve  detection  as  a  function  of  signal-ttMxnse  ratio  can  be  detennined.  Broadband 
detection  in  which  the  detection  statisdc  is  based  on  bispectnim  values  over  die  entire  principal 
domain  is  discussed  in  Hinich^*  as  opposed  to  narrowband  detectkm  at  a  single  point  in  the 
bispectrum. 

A  comparison  of  the  detection  performance  of  both  the  normalized  and  unnotmalized 
biqrectnrm  to  the  power  qrectrum  for  detecting  harmonics  in  Gaussian  noise  is  drown  in  the 
Rgure  2.0-1  fiom  Wilson.^  This  figure  is  based  on  tiiree  harmonics  havirtg  the  same  signal- 
to-noiae  ratio.  The  number  of  averages  used  was  10  necessitating  an  ^[rproximatimr  to  the 
a^mpiotic  statistics  used  to  produce  the  results  dnce  the  etqwessioo  for  die  exact  distribution 
has  not  been  dmived.  The  unnoimalized  bispectrum  (B)  perfotins  better  than  die  power 
plectrum  (PS)  idifie  the  normalized  biqiectium  (MB)  performs  worse.  The  best  that  can  be 
said  about  biqiectial  estimation  based  iqxn  this  ROC  (Receiver  Operating  Characteristic) 
curve  is  diat  it  is  no  better  dum  qiproximaiely  2*2.5  dB  over  power  spectral  techniques  in  the 
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signal-to-noise  ratio  range  of  -10  to  10  dB  with  a  normalized  version  up  to  8  dB  worse  on  the 
hi^  end.  It  is  worth  noting  that  Wilson  assumed  that  all  of  the  vaiiaUes  critical  to  successful 
bispectial  processing  are  optimal.  i.e..  sufiiciait  signal  skewness  exists,  the  noise  is  Gaussian 
or  ncm-Gaussian  symmetric,  and  the  processing  parameters  (DPT  length,  number  of  DFTs 
averaged,  etc.),  are  chosen  such  that  detection  will  result  from  a  sufficiently  large  noncottrality 
parameter  for  the  noncentral  chi-square  distribution  in  those  cases  where  the  skewness  function 
is  low. 


Figure  2.0-1  Comparison  of  Bispectrum  and  Power  Spectrum 

NaiTvwtMuid  Dtfection  Performance  for  PFA  s  OjOOI 

2.1  NONPARAMETRIC  METHODS 

The  methods  can  be  broken  down  imo  two  classes,  direct  and  indirect  The  methods 
have  the  advantage  of  being  strais^  forward  allowing  die  introduction  of  polyspectnd 
(Mocessing  fh»n  a  computational  standpoim  which  is  comparable  to  the  wdl  known  classical 
Discrete  Fburier  Transfbnn  (DPT)  implementaticm  of  the  power  spectrum. 


20 

In  practice,  conventional  bispectrum  estimators'^  use  a  finite  set  of  observatitm 
measurements.  However,  the  caution  is  that  limitations  on  the  gtarigtirai  variance  of  the 
estimates,  computer  time  and  memory  requirements  present  severe  imfdanentation  proUems. 

In  general,  the  direa  and  indirect  methods  for  topectrum  eaimation  give  differem  statistical 
estimates.  These  estimates,  however,  are  asymptotically  unbiased  and  consistent  and  have 
distributions  that  tend  to  be  comfdex  Gaussiaa^^  Conventional  estimators  have  hi^  variance 
and  thus  require  a  large  number  of  records  to  obtain  smooth  bispectral  estimates.  Increasing 
the  number  of  records  or  segments  causes  the  tradeoff  of  increased  computation  time  and  can 
present  additional  ixoblems  with  nonstationarities.  Frequency  domain  averaging  over  small 
rectangles  to  reduce  variance  has  the  unwanted  side  effect  of  increasing  laas  as  well  as 
computatian  time.  As  per  power  qrectrtim  computations,  the  effective  number  of  records  can 
be  increased  by  overiapping  for  short  data  records.^^ 

2.1.1  Direct  Method 

Hus  class  of  ocmventitHUd  estimator  is  useful  for  the  generation  of  rnomem  qiectra 
using  FFT  algorithms  and  is  defined  in  Nildas^^  as  follows. 

Let  {X(1),X(2),...  JC(N))  be  the  available  set  of  observatitHis  for  Uspectrum  estimaiiotL 
Assume  foat  f,  is  the  samfding  fiequency  and  »  fjff^  is  foe  required  spacing  between 

frequency  samples  in  foe  biqrectrum  domain  along  horizontal  or  vertical  directions.  Thus 
is  foe  total  number  of  frequency  samides. 

1.  Segmem  foe  data  into  K  segments  of  M  samples  each,  foat  is  N  &  KM,  and  subtraa 
foe  average  value  of  each  segment.  If  necessary,  add  zeros  at  each  segment  to  obtain 
a  oonveniem  length  M  for  the  FFT. 
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2.  Assuming  that  }  ue  the  dsu  of  segment  (i).  genenie  the 

IXT  ooefRdems 


jr-t 


«  s>« 


X  «  0,K..M/2 

i  = 


3. 


ai.M) 


In  general,  ii  •  x  where  »  a  positive  int^er  (assumed  odd  number), 
dut  is,  j/j  >  2£|  -t-  1.  Since  a  is  even  and  is  odd,  we  compromise  on  die 
value  of  (closest  integer).  Generate, 


E 

Vl--*! 


ai.1-2) 


fori  ■  1^.X  example,  in  die  qiedd  case  of  die  biqiectium  where  no 
averaging  is  perfonned,  dut  is,  Af|  *  1,  Lj  «  0«  ^  have  the  triple  produa 

A,*  a,.  X,)  -  -^1^(X,)  I^x,)  r»'{x,+x,)  c.i.1.3) 

^0 


Equadon  (2.1.1-2)  gives  the  user  the  qxion  of  peifonning  averaging  over  neighboring 
frequencies  to  reduce  the  In^Kctnim  estimation  variance.^  The  number  of  adjacent 
fie^iencies  over  which  diis  averaging  is  to  be  performed  is  determined  by  the 
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praneter  «2L|  *l  >nust  always  be  odd  so  ttiat  the  number  of  points  in  both 
diiectioas  and  parallel  to  the  biqxctnim  axes  will  have  a  mi(4)oint  or  median 
that  will  contain  the  averaged  value  from  both  directions  as  given  in  equation  (2.1.1- 
2).  Fm*  and  we  have  no  frequency  averagiqg  with  equation  (2.1. 1-3) 

resulting  and  *3  equation  (2.1.1-2)  would  reflea 

frequency  averaging  over  three  adjacent  frequency  points  along  both  axes  of  the 
btqrectnim  with  the  DPT  perfonned  over  points  for 

each  segment  Thus,  foe  desired  numba  of  frequency  samjdes  is  attained  as 

At  the  outsa  of  using  the  direa  mefood  we  specify  vfoich  is  fypically 
some  convenient  value  for  FFT  computation.  If  we  also  have  the  frequency  averaging 
option  available  to  us  then  we  can  iqxdfy  a  value  for  deflning  our  avoaging 

interval  lengfo  akmg  both  bispecoum  axes.  This  changes  our  original  value  of 

wliidi  now  becomes  which  is  our  new  or  effective  value  of 

reflecting  die  frequency  averaging,  Le.,  effective  <  original  Since  is  even 

and  ji  is  odd,  foe  new  value  of  is  usually  not  an  integer.  We  conqnomise  by 

dmosing  die  closest  integer  to  it  wMch  is  used  in  the  calculadon  of  the  scaling 
parameter  for  equation  (2.1.1-2).  The  averaging  opetadon  does  indeed  reduce  foe 
variance  but  it  may  also  introduce  Uas.^  This  brute  force  fretpiency  avoaging 
operadon  is  impfemented  in  ifi-Spec*s  direa  U^iectium  esdmadon  routine  "biqtec.d”^ 
via  a  2-D  oonvoludon  of  a  frequency  domain  smoothing  window  wifo  an  averaged 
bispectrum  over  if  segments  (records)  udddi  is  defined  next 
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4.  The  moment  spectrum  of  the  given  data  is  die  avenge  over  die  x  segmoits 


where 

An  optimal  frequency  domain  smoothing  window  in  die  mean  square  enor  sense  (MSB),  the 
Rao-Oahr  window*^  which  optimizes  the  tndeoff  between  variance  reducticn  and  bias 
introduction,  is  offered  as  a  Hi-Spec  option.  The  Hi-Spec  window  imidementation  is  defined 
as  fidlows:^ 

1.  The  panmeier  wind  specifies  die  fiequency-domain  smoothing 

window.  If  wind  is  a  scalar,  the  Rao-Gabr  window  of  lengdi  wind  will 
be  used.  This  window  is  defined  by.*^ 


W(,m,n) 


*  *  mn 

N* 


(ir,r)cG 


where  N  is  half  the  FFT  lengdi.  iffft,  and  G  is  the  set  of  ptdnts,  (m,  a), 
smisfying  the  ellipse 


m 


2 


+  H* 


mn<  *■ 


wmS 


2.  A  unity  value  fiir  wind  results  in  no  windowing. 

3.  If  wind  <s  0.  die  defuilt  value  of  S  will  be  used. 
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4.  If  wind  is  a  vector,  it  is  assumed  to  specify  a  1-D  window  fitom  which 
a  2-D  window  is  computed,  W(mjt)  •  w(m)w(H)w(m*H). 

5.  If  wind  is  a  2-D  matrix,  it  is  assumed  to  specify  the  2-D  smoother  diiectiy. 
The  net  effect  of  this  windowing  procedure  for  stqis  1  and  3  is  that  all  of  the  ptrints  G.  the 
shaded  hexa^mal  region  in  Figure  2.1.1-1  fmn  Rao^^  comprising  the  12  symmetry  regions  of 
the  biqrectnim  (see  Appendix  B)  bounded  by  the  ellipse  area  Cp  are  windowed  where  Gj  is 
die  dlipae  area 


G  is  the  rinded  hexagonal  area 


+  aw  <■ 


wind* 


(a^R):  N  N  *  <»2T"jr^T» 

VVVZ) 

and  G  <  G/. 

The  optimal  Rao-Gabr  window  reflects  the  minimization  of  the  MSE  consisting  of  the  variance 
of  tfK  bispectial  estimate  plus  the  square  of  its  bias.  This  weight  function  is  actually  an 
etqxesskn  contained  in  the  variance  of  the  estimate  and  is  riiown  to  be  smaller  than  ar^  other 
biq)ecttal  weight  function,  dieteby  minimizing  the 
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figure  Rao^abr  Window  Bounds 

2.14  Jl  Test  Results 

As  described  in  Appendix  C  synthetic  signals  at  vaiying  signal-io-noise  ratios  were 
generated  to  exercise  the  direct  biqiectnim  algoridun.  Results  for  synttieiic  signals  embedded 
in  Gaussian  white  noise  indicate  no  significant  impiovemeitt  in  biqrectial  analysis  over  power 
spectral  analysis.  The  reason  fbr  diis  is  that  the  signals  simulated  exhibited  a  symmetric  pdf 
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(riBewaeii  inmn.  0)  wliidi  made  it  difBcult  for  die  aigoridan  to  difieiendate  between  non- 
Cfonmian  qnoBmetric  sigii^  nd  non-Gaumiaa  qmunet^  Bodi  of  theae  are  equally 
auppreaaed  aa  the  baalc  miimptinna  in  regard  to  biqiecsial  analyaia  are  that  the  signals  are 
noo>Oaumian.  zeio-men  with  non-zeio  skewness.  i.e..  exhibit  non-symmetiy.  «d  thm  the 
noise  is  KiD-mean  Gaussian  or  non-Gaussian  symmetric. 

Real  acdve  sonar  field  data  was  processed  with  similar  results  to  the  simulaied  data  for 
similar  reasons,  U..  die  projected  signals  used  for  pingirig  exhibit  symmetric  pdft  udiich  were 
unperturbed  upon  backscatlering  fiom  a  target. 

In  sununary.  the  results  indicaie  that  fai^iectral  techniques  which  siqipress  Gaussian  (or 
nop-Oaussian  •  symmetric  pdf)  noise  suffer  fiom  a  lack  of  signal  skewness  which  precludes 
any  discemabie  huprovement  over  power  detection  methods  (see  section  C.2  of  Appendix  C). 

2.U  fiadbed  Method 

TUs  dass  of  oonvcndonal  estimator  is  based  on  the  2-D  ITT  of  a  windowed  third 
order  cnmnlanc  fimedon  and  is  defined  in  Nikias**  as  fidlows. 

^  be  fibren  dam  set  Then  we  have  the  following. 

1.  S^ment  the  data  inro  g  records  of  j#  samples  eadi.  that  is,  /f  •  ru. 

2.  Subtrad  the  average  value  of  each  reomd. 

3.  Assuming  that  -  OlU.Jir-1)  is  the  dam  set  per  segmem  j  .  i, 

2^.,  g,  obtain  estimates  of  the  Mgher-order  moments 

4.  Average  over  an  segounts 

5.  Generare  die  irdHirder  cumulant  sequence  b  is  a  function  of 


1  ^ 

Xf  M, 

n  ■  23..^  <  ■  (2.12-1) 

Sj  ■  10&X(0(~T  !•••••  _|), 

^  -  ndiiCJf  -  1,1#  -  1  -  T„J#  -  1  -  -  1  -  T,.,), 

I  X 

^•(Xif~»%.|)  *  "p  ^  (X|»...#t^.j)  (2.12-2) 

Ji  »  23.~**«  -  1. 

mamems  fimn  Older  fecood  to  nth  (see  RoeenNstt*  ^  Ibr  geneial  itlationsliip). 
However,  fiir  a  zero-fnean  process  we  have 


and 


•  I^(X|»T2,T3)  “  43(‘^iV^(X3  “  tj) 
-  li^*(T,X(T,  -  T^. 


(2.12-3) 


6.  Genentfe  the  Uglier-Older  ^Kctun  estimate 

L  L 

- 1)  *  ♦-  ^/(^l*****^*-!) 

»,  -  -L  *,  -  -t 


ai.2-4) 


xibtit  LkM  -  1  and  .  t)  ^  *  window  function.  Instead  of  a  cumuhnt 

spectrum,  a  moment  spectrum  can  be  genenued  1^  equation  a  12-4),  if  we  substitute  the 


estimates  given  by  equation  (2.12-2). 


Hie  multidimenskMial  window  function  ifaould  otisfy  the  symmetry  piaiwfties  of 
M^to’-oider  moments  or  cumulants.  They  should  also  be  zero  outside  the  of 
siypMt  of  estimated  cumulants  and  have  noonegative  Fouiig  transfonns.  Adassof 
functions  that  satisfies  these  constnints  is  the  ftdlowing. 


udieie 

dCx)  -  dC-t) 

d(T)  ■  OlT>I. 

dm  «  1  (2.U-5) 

L 

d(T)esp(->mT)^  0,  fiir  aU  ta. 

Identity  (2.1^>S)  allows  a  reoonstniction  of  window  functions  for  faigher-oider  specoum 
estimation  using  standard  one-<limensional  (1  -  4)  windows.  However,  not  dl  i  -  4 

windows  satisfy  constraint  D(o)2  0  foe  aD  m.  Two  windows  witti  good  performance  in 
terms  of  bias  and  variance  minimization  are  (see  Nikias^^  and  references  tiierein): 


%  L  L  L 

0. 


ItI^I 

|t|>L. 


(2.1.2-6) 


and 
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I  <|tUt 


au-7) 


la  |t|>i. 

The  fifit  window  actdeves  a  bias  that  is  ^>oat  IS  peiooit  smaller  than  that  of  the  second 
(Panen)  window.  However,  the  fint  window  produces  variance  that  is  about  26  percent  laiser 
than  thst  of  Panen  window.  Equation  (2.1J2-7)  is  applied  in  Hi-^tec  to  minimize  the 


variance. 


2.1  Teat  Results 

Conventional  indirect  biqpectnl  analysis  employing  2-D  FFTs  of  higher  older 
cumulants  were  examined  lor  conqMiison  with  direct  methods  with  the  result  that  resolinion 

was  fitf  superior  using  the  direct  algoritlm  (see  section  C.1  of  Appendix  C).  The  inability  of 

( 

the  indirea  biHXCtnim  lo:  (1)  resolve  quadratic  phase  cotqded  components  rimOariy  to  the 
direct  bispectrum;  and,  (2)  suppress  uncotqded  sinusoids  mixed  widi  coiqtled  sinusoids;  caused 
questions  to  be  laiaed  as  to  its  ^ecthreness  Ux  noise  suppression  (uncoupled  sinusoids  should 
look  like  non-Oaussian  symmetric  pdf  noise  to  the  biqiectnim).  Thus,  indirea  biqwctia 
comparison  to  power  spectra  was  dropped  pending  a  mme  detailed  examination  of  die  dieoiy 
and  implementation  of  conventional  indirect  bispectrum  methods.  (See  section  C4 
DIRECT/INDIRECT  RISPECTRUM  COMPARISON  FOR  SIMULATED  DATA.) 
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Z2  PARAMETRIC  METHODS 

An  altemative  appioadi  to  Fourier  based  methods  for  the  nueipretation  of  time  series 
data  is  the  oonstnictioa  of  white  noise  driven  linear  parametric  models  from  the  undolying 
physical  process.  The  motivations  fior  this  are  thieefirid:  (1)  to  recover  phase  infonnation 
accurately,  (2)  to  increase  the  lesdution  capability  of  an  estimator  involving  closely  spaced 
peaks  in  the  bispectium  domain,  and  (3)  to  increase  bispectrum  fidelity  in  the  sinuukms  where 
non-Oaussiao  processes  are  actually  parametric  or  may  be  approximated  by  parametric 
models.*^  These  mediods  typically  suffer  in  low  signal>«M90ise  ratio  environmems  and 
demand  a  certain  amoum  of  qxiori  knowledge  tiriudi  isn't  always  readily  availaUe  in  practice. 

Parametric  methods  firr  U^ier-order  spectrum  estimation  are  described  in  section  7.4 
of  hnidas.*^  MA.  noncausal  /\R.  and  ARMA  methods  are  discussed  reqrectively  in  relation  to 
a  nonminimum-phase  system  identification  problem  ^di  staled  simi^y  is,  given  only  a  finite 
length  of  data,  it  is  required  to  estimate  the  bi^wctrum  of  the  underlying  discrete  random 
process  via  parametric  modeling  of  its  third  moments.^^  For  example,  the  MA  method 
requires  the  estimation  of  tiie  coefificiems  of  the  MA  process  via  a  epical  second-order  metiiod 
to  reflea  die  autocorrelation  structure  of  the  data.  This  apriori  knowledge  is  then  ap{died  to 
the  estimation  process  with  die  end  result  being  a  fine  tuniitg  to  the  coarse  estimate  made  via 
die  second-order  process  to  improve  resolution.  A  typical  example  of  this  in  Hi-Spec  is  the 
"qpc.lM”  program  defined  under  quadratic  phase  coiqiling  in  section  A.0  of  Appendix  A.^ 
This  program  uses  a  "third-order  recursion”  to  estimate  the  bispectrum  as  a  function  of  the  AR 
model  order  which  should  be  greater  than  the  number  of  qieoral  lines.  In  particular,  closely 
qMced  HKCtral  lines.  The  motivation  ftir  such  parametric  teOmiqiies  is  to  (1)  more  accuratdy 
estimate  biqiectra  of  an  AR  process,  nd  (2)  rmolve  dosely  qnced  peaks  in  die  biqiectrum  or 
delea  die  presence  or  absence  of  quadratic  ptux  coiqiling  at  frequency  pdrs  in  dose 
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pRHdmity  to  each  othor.  The  conjecture  is  made  that,  just  as  AR.  ARMA  models  offer  high 
resolution  capabilities  as  altenutives  to  conventional  power  q)ectium  methods,  dieir  tMq)ectial 
counterpaits  can  be  expected  to  offer  similar  hi^  resolution  capaMlities  ctunpared  to 
OMiventional  direa  or  indirect  bispectral  estimators.  Comparison  examines  are  given  in 
Raghuveer.^ 

There  ^)pean  to  be  some  difference  of  (^xnioo  in  the  literature  conceining  die  efficacy 
of  assuming  that  a  given  time  aeries  can  be  described  by  a  low  order  AR  or  ARMA  model  as 
iqmited  in  Thomson^  which  details  die  authois*  eiqierience  in  analyzing  many  differem  time 
series  models  from  sdetdific  and  engineering  pitdilems  with  die  statement  that  "neidier  of  us 
has  yet  seen  an  example  iiriiich  could  leasonatdy  be  described  by  an  autoregressive  model" 
This  statemett  is  later  qualified  by  die  admission  duo.  while  stating  dut  such  examples  may 
exist,  the  probability  of  realizing  one  appears  to  be  so  low  diat  "starting  widi  the  assunqition 
that  a  given  time  series  can  be  described  by  a  low-order  AR  or  ARMA  model  is  a  prescription 
for  trouble."  This  is  not  the  teascm  for  neglecting  parametric  biqiectral  esdmators  in  this  paper 
as  mudi  as  the  requirement  for  the  minimizadtm  of  tqnimi  knowledge  in  regard  to  bispectral 
signal  detecdon  versus  sectmd-order  stadsdcs  power  spectrum  techniques,  ix..  the  main  focus 
here  is  coarse  as  opposed  to  fine  tuning  altbou^  some  altemadve  methods  emidoying  hi^- 
order  stadsdcs  for  resoludon  improvement  will  be  listed  in  section  22.1. 

2JL1  Harmonic  Retrieval  Test  Results 

As  an  adjuiKt  to  the  problems  encouraered  in  Nspectral  signal  detection  due  to  the 
lack  of  edio  or  radiated  signal  skewness  compared  to  noise  ctmtributitms  having  negligiUe 
skewness,  frequency  resoludon  tests  posatde  using  die  Hi-Spec  routine  "harm.est"  woe 
performed  usirtg  signals  generated  by  "harm_gen"  (see  Appendix  A).  The  "harm.est"  routine 
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gives  an  estimarion  of  the  frequencies  of  real  hannonics  in  noise  witfi  comparison  to  a 
convettional  periodogram.  The  estimation  is  given  using  the  MUSIC,  Eigenvector.  Pisareidu), 
ML  (Capon)  and  AR  meduds  based  on  the  diagonal  slice  of  fourtlKMder  cumulants.  No 
improvement  was  noted  over  nonparametric  FFT-hased  direa  biiqKctrd  algmiduns  as  a  FFT 
biqiectral  slice  exhibited  the  same  peak  to  peak  reqxmse  as  the  diagonal  fourth-order  cumulant 
slice.  I^Meover.  for  low  signal-to-noise  ratio  cases,  these  higher-oider  statistical  resolution 
improvement  methods  deteriorated  significantly  in  comparison  to  the  classic  periodogran  -  a 
problem  common  to  their  second  order  coumerparts. 

It  is  noted  diat  coarse  adjustmem  Omg  range  detection)  as  opposed  to  fine  (target 
classification)  is  the  main  focus  of  this  paper,  b  was  the  audior*s  iruent  to  call  attention  to  the 
fact  dtat  Ugher-order  methods  are  no  resolution  panacea.  Lack  of  signal  skewness  inhibits  any 
propagation  gain  as  die  signal  is  suppressed  etpially  with  the  noise.  Lack  of  sufiiciem  SNR 
inhibits  any  restdudon  improvemenL  The  latter  is  not  a  new  resub.  The  fbnner  verifies  the 
examples  fiom  die  literature  review  and  the  tests  performed  for  this  paper  as  rqwrted  in 
^ipendixC. 
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CHAPTERS 

CONCLUSIONS/FUTURE  WORK 

This  diapier  examines  the  efficacy  of  signal  detection  in  underwater  acoustics  using 
polyspectral  mediods  in  the  light  of  current  research  into  nonstationary  HOS  with  future 
directions  for  study  indicated. 

GENERAL  CONCLUSION 

The  theme  throughou  this  paper  is  that  the  critical  panuneter  in  relation  to  the  third 
Older  cumulaitt  spectrum  (bispectntm)  as  defined  in  equation  (B.2-S)  successfully  realizing  any 
type  of  propagation  gain  over  the  second  order  cumulant  spectrum  (power  qrectrum)  as 
defined  in  equation  (B  J*3)  is  the  signal  skewness  as  defined  in  equation  (B.1-6).  Prelimmary 
results,  pending  a  more  detailed  examirution  of  simulatedAeal  data,  indicate  that  if  the  echoed 
(active)  or  radiated  Qrassive)  signal  does  not  exhibit  a  degree  of  moderate  skewness  which 
remains  to  be  quantified,  but  as  a  rule  of  thumb  approximately  .5  where  1  indicates  a  hi^y 
skewed  distribution,  the  probability  that  HOS  agnal  jHocessirig  sr^rhistication  will  realize  a 
sigidficant  d^ection  improvement  over  sectmd  order  methods  iqrpears  to  be  low.  This  is  borne 
out  1^  observing  that  a  time  aeries  representing  a  mixture  of  quadratically  phase  cou{ded 
sinuaoids  and  uncoupled  sinustdds  exhibits  a  modentte  skewness  measure  while  a  similar  series 
iqriesenting  only  unoouided  siraisoids  does  not,  i.e.,  suppression  of  the  uncoupled  siiuisoids  is 
realized  in  die  former  case  virile  there  is  no  difference  between  uncoiqded  sinusoids  and 
Gaussian  or  nonOaussian  symmetric  noise  in  the  latter  case. 
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3.1  NONSTATIONARY  HIGHER-ORDER  SPECTRA 

Success  in  the  area  of  nonstttionaiy  higher  order  spectral  estimators  resulting  in  a 
claimed  3  dB  detection  perfonnance  improvement  under  realistic  ctmditions  M^ien  cmnpared  to 
the  stationary  power  q)ectrum  has  been  reponed  in  Wilson.^  In  this  report  die  detection 
performance  of  different  types  of  higher  order  qiectra  are  discussed  for  a  vari^  of  signals. 
The  introduction  of  what  the  authors  refer  to  as  a  "new  type  of  higher  order  spectra  called 
nonstationaty  higher  order  ^xctra"  is  made  with  the  distinction  that  these  nonstatiotuuy  higher 
order  qwctra  are  not  stationary  higher  order  spectra  representations  of  nonstationaty  processes 
but  rather  are  different  qwctra  which  omtain  the  stationary  higher  order  spectra  as  a  subset  of 
thdr  domain.  It  is  shown  quantitatively  throu^  theoretical  predictions  and  simulations  that 
these  qiectra  perform  better  at  detecting  nonstationaty  agnals  tiian  do  die  traditional  stationary 
spectra.  The  detection  perfonnance  of  die  bispectrum  is  examined  in  comparison  to  the  power 
qiectrum.  The  motivation  for  the  introduction  of  this  new  class  of  qiectra  is  the  detection  of 
nonstationaty  signals  such  as  hamionically  related  narrowband  componetas  of  finite  duration 
transients.  A  fimctional  approach  is  used^  to  define  the  stationary  order  spKXn  or 
poly^Kctra  and  a  statistical  tqiptoach  is  used  to  define  the  nonstationaty  hi^r  order  spectra 
or  cumulant  spectra.  The  estimation  of  higher  order  spectra  are  discussed  firom  both  stationary 
and  nonstationaty  perspectives. 

The  stationary  fiinctimial  tqiptoach  defines  a  i»ocess  with  orthogonal  increments  based 
rm  die  classical  derivation  as  defined  on  page  32  of  Priestly.^  The  nrmstationaty  statistical 
approach  assumes  a  vector  valued  discrete  time  series  whose  ndi  order  cumulant  exists  and  is 
finite.  Then  the  nth  order  discrete  Fourier  transform  of  die  cumulant  exists  and  an  ndi  ordm* 
cumulant  spectral  density  function  is  subsequently  defined  as  the  mh  order  discrete  Fourier 
transform  of  .the  cumulant  function.  The  key  here  is  diat  no  assumption  of  statirmarity  has 
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been  made.  For  the  estimation  of  nonstationary  higher  order  spectra  die  result  that  the 
cumulant  of  die  finite  Fourier  tiansfoims  is  propoitional  to  the  nth  order  polyspectial  denrity 
fimcdon  |dus  lower  order  terms  is  extended  to  being  propoitional  to  the  mh  order  cumulant 
spectral  density  function  (defined  above)  lower  order  tenns.  The  observation  is  made 
that,  for  stationary  processes,  the  cumulant  spectnim  is  zero  except  in  the  domain  for  which 
the  sum  of  all  fiequency  compments  are  zero.  Thus,  within  this  domain  the  cumulant 
qwctnim  is  equal  to  the  dasacal  ptriyspectnim  defined  by  the  stationary  functional  ^ipioach. 
(This  domain  is  the  dark  shaded  triangular  region  in  die  block  labeled  "I”  in  the  first  quadrant 
of  Figure  BJ2-2  entitled.  Fundamental  Domain  of  the  Discrete  Time  Bispcctrum,  from 
Section  B2.  The  light  shaded  triangular  region  adjacem  to  this  domain  becomes  a  region 
where  non-zoo  values  indicate  nonstationarity.) 

3.U  Condusion 

The  Mqpectial  results  have  previously  been  discussed  in  this  paper  (see  Section  2.0). 
The  message  of  the  statistical  approach  taken  in  Wilson^  is  dut  Mindly  computing 
polyqiectFB  is  folly  widxwt  taking  the  time  to  understand  the  underlying  statistics  of  the 
processes  involved  triuch  necessitates  a  certain  degree  of  statistical  prqxocessing  prior  to 
attempting  any  type  of  polyspectial  analysis.  The  nonstationary  statistical  ^iproadi  used  in  the 
rqmit  to  define  cumulant  spectral  density  functions  lends  itself  to  a  better  understanding  of  the 
statistics  involved  in  pMyspectial  {Hocessing  and  is  an  example  of  the  multidimensional 
statistical  basis  used  for  tri^iectrum  (kurtosis)  conputations  diat  will  be  addressed  in  sectioi 


3.2. 
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3.U  Fkiturc  Work 

Fuittier  examinatian  is  lecommended  to  better  understand  the  nuances  between 
stadonaiy  and  nonstationaiy  higher  older  spectra  computatioas  ftom  both  algorithmic  mi 
interpretive  viewpoints  as  a  function  of  the  underlying  statistics  (real  and  assunwd)  of  the 
processes  involved. 

3^  TRISPECTRUM 

A  natural  extension  to  the  bispectrum  is  the  examination  of  fourth-order 
(kurtosis)  in  reference  to  determining  the  degree  of  success  for  signal  detection  in  die 
tiansfitMm  (trispectral)  domains.^^  Work  is  specifically  being  directed  toward  computing  the 
statistics  for  the  trispectrum  similar  to  that  for  the  bispectrum^*^*  over  tiie  region  of  support  in 
tile  first  octant  for  trispectral  computttions.  ix..  the  principal  domain  of  the  trispectrum  for 
stationary  continuous  time  processes  which  is  a  triangular  cone  in  the  positive  octant— the 
three  sides  of  titis  cone  being  die  itianes  fmmed  by  the  intersections  of  three  of  die 
tiispectium's  symmetry  fdanes.^^ 

As  defined  in  (B2-8)  the  triqiectrum  is  the  fourdi  order  cumulant  qrectrum  of  a 
random  process.  The  motivation  for  investigating  the  triqiectrum  as  described  further  in  Dalle 
MoU^  is  diat  bispectral  tests  don’t  serve  as  cmnidete  tests  for  the  rgectkm  of  Gaussianity 
and  linearity  hypodieses.  Althoush  impractical,  a  conqdete  lest  based  on  hi^ier  order  qiectra 
would  involve  all  possible  prdyqiectia  as  non-Gaussian  processes  may  exist  vriiich  have  zero 
biqiectial  values  and  in  turn  zero  dtewness  measure  over  the  comidete  principal  domain.  It  is 
observed  that  since  most  non-Gaussian  processes  have  non-zero  hi^  order  cumulants,  it  is 
unlikdy  diat  a  process  would  have  both  its  skewness  and  kurtosis  functitms  identfoally  equal  to 
zero  over  the  entire  principal  domain.  Therefore,  just  as  additiorud  statistical  goodness-of-fit 
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ten  preeent  stronger  cases  Im-  aoceptanoe/rajeoion  of  hypotbeses.  bispectial  and  trispectial 
baaed  tests  togetber  make  a  stronger  case  for  die  lejecdon  of  the  null  bypothesis  of 
Gaussianity.  Similarly,  tbe  constancy  of  the  skewness  function  is  not  always  a  good  indicator 
for  a  linear  model.  Thus,  die  tiiqiectrum  presents  a  more  powerful  tool  for  the  nalysis  of 
non-linear  vid  non-Gaussian  time  aeries.  For  the  tri^KCtrum.  the  square  of  the  avenge 
kurtosis  funr  r »  becomes  the  test  for  Gaussianity.^  The  symmetries  vdiidi  need  to  be 
considered  in  order  to  define  the  continuous  and  discrete  time  principal  domains  of  the 
trispectrum  along  widi  die  deaciiptian  of  die  wedge-shaped  hyperplane  cone  m  the  frequency 
triple  (flXUS)  which  forms  die  continuous  time  principal  domain  are  given  in  Dalle  Mtdle.^ 
This  pyramid  shaped  wedge  is  defined  exjdicidy  Ux  band  limited  processes.  Tbe  principal 
domain  of  the  trispectrum  for  a  discrete  time  band  limited  process  is  derived  in  a  similar 
manner  to  die  bij^KCtrum’s  principal  domain.  Samplmg  at  the  Nyquist  frequency  results  in  the 
discrete  time  principal  domain  being  a  pyramid  widi  a  triangular  base  that  is  larger  than  the 
support  set  for  the  ooittiiuous  time  band  limited  trispectrum.^ 

Kurtosis  esdmation  has  also  been  extmsively  addressed.^*^  The  occurrence  (tf  non- 
Gaussian  dgnals  in  underwater  acoustics  doe  to  mobipadi  and  modulation  effects  is  the  main 
motivation  here  as  sinusoidal  and  narrowband  Gaussian  signals  idddi,  tdwn  propagated 
dirou^  feding  w  muldpmh  environments,  are  received  as  non-Gaussian  in  terms  of  the 
frequency  domain  kurtosis  estimate.^  The  results  are  rqjfdicable  to  bodi  active  and  passive 
sonar  with  die  active  case  conjectured  to  be  non-Gaussian  (dcewness  at  the  returns)  due  to 
modulation  effects.  Results  have  sug^sted  the  possibility  of  detecting  non-Gaussian  signals 
via  kurtosis  estimation  at  lower  signal-to-noise  ratios  than  the  power  plectrum  method. 

For  under  ice  data  the  non-Gaussian  signals  are  due  to  ice  movement  which  produced 
transient  and  modulated  signals  in  die  passive  cme.^  (If,  in  feet,  suffidem  skewness  is 
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obierved  the  retums  for  the  active  caae  then  third  oida'  stadstici  should  geneiaDy 
diaaiminaie  between  non-Gaussisn  and  Gaussian  signals.^  Kuitosis  is  defined  to  be  die 
ratio  of  a  fourth  order  central  moment  to  the  square  of  a  second  order  central  moment  It  has 
been  shown  to  be  a  locally  optimum  detection  statistic  under  certain  conditions  and  is  extoided 
in  the  frequency  to  the  detection  of  outlien  horn  an  otherwise  Gaunian  sample  -  the 
outliers  being  equivalent  to  the  randomly  occurring  signal  that  it  to  be  detected. 

Based  on  the  analyses  of  real  underwater  acoustics  data,  conditions  exist  under  whidi 
frequency  domain  kuriosis  estimation  indicates  die  presence  of  randomly  occurring  signals.^^ 
For  a  Gaussisn  distribution,  kuitosis  will  have  a  value  of  around  3  widun  some  bound  as  a 
function  of  the  number  of  samples.  For  randomly  occurring  signals  dut  produce  non-Gaussian 
distributions,  die  kurtods  esdmaie  can  be  less  than  or  much  greater  than  3.  A  modd  for  the 
received  data  which  contained  die  effocts  of  amplimde  and  phase  fluctuation  of  the  signal 
along  with  modulation  effocts  was  used  to  evaluate  die  kuitosis  estimates  in  Dwyer.^^  The 
ftequency  domain  kurtosis  estimadon  was  computed  and  diown  to  have  significantly 
values  even  for  small  agnal-to-noise  ratio  cases  surpassing  the  averaged  power  qiectium 
estimation  in  some  instances  due  to  its  smsidvity  to  the  pdf  of  die  signal  Results  flom 
Dwyei^  suggest  that  for  randomly  ocoming  signals.  e.g..  a  transient  or  a  modulated  signal 
over  long  integration  times,  die  ftequency  domain  kuitosis  estimate  may  be  a  better  detection 
stadsdc  diat  the  power  spectium  esdmaie.  For  fading  environments  wluch  frequendy  occur  in 
underwater  aoousdcs.  die  kuitosis  estimates  are  significandy  eriumced. 

The  udli^  of  die  fourth  order  qiectium  for  the  dassiflcadon  of  transient  signals  in 
passive  sonar  and  resonances  in  active  sonar  along  with  the  aunpotadoiial  drawbacks  is  given 
in  Dwyer.^^  The  importance  of  diis  is  obvious  because  many  noise  sources.  i.e..  ambient, 
revetberadon.  flow,  etc.,  are  Gaussian  compared  to  signals,  i.e.,  dnusoids,  transients,  active 
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miff  tnmmifiiont.  wlridi  are  all  non-Oaussiaii.  It  is  demonstrated  that  the  fbuith  mder 
spectrum  diffetenliaies' between  a  sum  of  sinusoids  and  duis  can  aid  in  active  sonar 
classificatioa  When  long  pulse  trains  are  transmitted  the  return  echo  is  usually  modulated 
fram  taiget  dynamics,  Doppler  spreading,  or  medium  effects.  Under  these  conditions  the 
plectrum  can  be  severely  distotted  making  detection  and  dassification  impossible.  Using  an 
coded  pulse  train  it  is  shown  dutt  die  fouith  order  spectnim  can  extract  range  and 
Doppler  taifonnadon  wUk  for  die  same  condidons  the  spectnim  is  uaeless.^^  The  emphasis  is 
on  the  special  case  of  die  ftmith  order  cumulam  defined  as  kuitosis  in  (B.1-6).  The  treatment 
here  is  multidimensional  as  a  new  class  of  non-Gaussian  density  functions  which  iqxesent 
mcanjngftil  signals  Qndependent  data  are  not  required)  is  introduced  with  their  conesponding 
fourth  order  spectra.  The  importance  of  this  work  is  idated  to  the  stadsdcal  treannent  of 
modulated  processes  as  related  to  sonar,  in  particular  active  sonff . 

To  demonstrffe  theoretical  results  in  r^ard  to  frequency  domain  kurtosis  estimation  an 
experiment  was  conducted^  vdiere  a  sinnsoid  was  modulated  by  white  Oausaian  noise, 
transmitted  dirough  die  water  and  received  on  an  omni-directional  hydrophone.  TUswasan 
attempt  ff  simulating  the  modulation  due  to  target  effects  witich  was  previmidy  discussed. 

The  second  mder  spectnim  (power  spectnim),  the  qiectium  of  the  special  case  tff  die  fouith 
order  moment,  and  the  qiectrum  of  the  spedal  case  of  the  fourth  order  cumulsnt  (tiiapectrum) 
derived  at  the  omput  of  a  lowpass  filler  were  estimated  from  die  filtered  data.  The  experiment 
verified  dieoretical  results  by  showing  that  die  second  order  qxctrum  could  not  extract  the 
sinusoidal  fnequency,  but  the  reflective  fiectrums  of  the  fourdi  order  moment  and  cumulant 
oouhL  It  is  further  reported^  dut  the  analysis  of  sonar  data  clearly  diow  modulated  or 
fluctuating  signals.  As  the  range  increases  between  a  source  and  receiver  diese  fluctuating 
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of  the  tnmmitted  signal  in  tbe  experiment 

The  piopefties  of  the  special  case  of  a  fouitfa  older  cumulant  m  oonsideied  via 
impkmematioa  using  1-D  Fourier  transfiMms  as  opposed  to  the  computational  complexity 
invtdved  with  ftiU  3-D  Fourier  transfonn  realizations  of  the  trispectnim  as  defined  in  (B  J-8). 
The  derivation  of  the  Fourier  transft»m  of  the  special  case  of  a  fiMiith  older  cumulant  ft>r  the 
extraction  of  range  and  Doppler  inlbimation  is  given  in  Dwyer.^^*^^  An  improvement  on  the 
Older  of  that  noted  in  Dwyei^  for  die  detection  of  siniiaoids  was  similarty  observed  when 
oonqiaring  fourth  order  cumulant  lechnkpies  to  the  second  order  spectrum  for  range/Doppler 
Implications,  ft  was  demonstrated  by  simulations  diat  the  Fourier  tmsfoim  of  a  qiecial  case 
of  die  fourdi  oider  cumulant  could  extract  range  and  Doppler  information  with  high  resolution 
even  when  Gaussian  noiae  was  added  lo  the  edio  while,  under  the  same  condidoiis,  the 
ambiguity  ftmction  and  second  order  qiecttum  processing  showed  significant  Doppler 

Conduaioa 

It  is  recognized  dut  the  additional  computational,  statistical,  and  physical  meaning 
complexity  concomitant  widi  tri^iectral  compotations  may  be  the  deciding  factors  in  any  kind 
of  a  tradet^  with  a  confectuitd  improvement  over  ptdyqtectral  methods  of  lower  order. 


3JJ  Future  Work 

If  special  cases  of  fourth  order  cumulant  spectrum  (trispectium)  computations  can  be 
realized  which  mhiimize  the  cooqadatiooal  effort  required  for  odierwise  fiiU  triqiectrum 
realizations  dien  the  efforts  dted  above  merit  fiirdier  scnuiny  in  regard  to.  at  die  very  least. 
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signal  daisificaiiao  wtiteh  exploits  tlie  phase  mienctioas  of  higher  order  aariair*,  ni  agnal 
detection  per  se  in  comparison  to  second  older  mediods.  Abbreviated  kuitosis  estanatkm 
algofitfams  might  even  present  viable  attemadves  to  tfaiid  order  mediods  when  signal  skewness 
is  not  a  problem.  sufficient  skewness  exists  to  employ  bispectnl  tedmiques.  This  did  not 
appear  to  be  the  case,  however  for  the  red  test  data  examined  in  this  repoit  as  signal 
drewness,  specifically  lack  of  same,  presented  a  very  real  problem  in  r^aid  to  the  assumptions 
neoessaiy  for  suocessftil  applicatian  of  higher  mder  qxctia  for  detection  purposes. 

Dwyer’s  woifc  points  out  die  foct  that  power  spectrum  density  (PSD)  estimation  is 
essentially  a  second  order  measure  udddi  is  not  sensidve  to  the  statistical  namre  of  the 
signals.^^*^  He  offen.  as  an  ahemative.  simple  fourth  order  movement  (kuitosis)  teats  in 
both  fieqpiency  and  time  as  supplements  to  PSD  estimates  and.  in  some  cases  of  practical 
iinponanoe.  as  superior  detection  statistics  to  PSD  tedudques.  It  is  hiflhly  recommended  diat 
the  success  whidi  dwyer  adiieved  using  fieqoency  domain  kurtosis  methods  be  extended  to 
the  time  domain,  as  he  suggests,  for  both  active  and  passive  sonv  scenarios.  Dwyer 
emphasizes  that  the  time  domain  is  nothing  more  than  a  special  case  of  the  frequency  domain 
and  diat  nalogous  results  should  also  hedd  in  the  spmial  domain  where  kurtosis  methotto  could 
be  evaluated  in  r^ard  to  target  angfc  estimation.  Ifis  past  success,  for  low  input  signal-to- 
noise  ratio  cases  (SNRi)  using  kurtosis  methods  exclusive  trf  any  cumulant  connotations  where 
kurtosis  methods  detected  non^Saussian  signals  tddle  PSD  methods  foiled  to  do  so.  is  reason 
enough  to  verily  his  conjectures  in  both  time  and  spatial  domains. 

Following  Dwyer’s  lead,  algofittuns  need  to  be  developed  t^ndi  examine  both  the  real 
and  hnaginaty  parts  of  foe  input  data  and  optimaDy  select  die  part  or  post  processing  of  both 
parts  wUdi  gives  die  best  detection  statistic  as  a  function  of:  (1)  die  signal  set  used.  (2)  die 
target  aspect  angle,  and  (3)  die  kuitosis  window  kngth  which  defines  a  diding  kuitosis 
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window  aciMS  n  acoustic  cyck  (ping  in  acdve  cate.  liAeniiig  interval  in  passive  The 
alforittams  need  to  be  exsmiiied  in  legaid  to:  (1)  edge  effectt  lewldng  fiom  the  windowing 
pioceduie  ft>r  simulated  high  and  low  kunosis  signals,  ^tedfically,  how  edge  effects  can 
maximize  detectkm  by  dictaiing  whether  tdgh  or  low  butotis  is  sought  for  the  detection 
statistic.  (2)  the  effect  ot  inducing  high  or  low  kuitosis  to  a  signal  distribution  prior  to  pinging 
for  the  active  case,  specifically,  how  sudi  inducement  effects  the  signal  upon  backscattering 
fiom  a  taiget  with  particular  attention  paid  to  the  degree  of  kuitosis  observed  post  target 
dynamics  andfor  modulation  effects;  (3)  evaluation  against  •««nuiat»d  test  beds  coittaming 
realistic  targets  from  widdy  acoqtied  sonm  models  as  a  function  of  taiget  aspect,  platform  and 
taiget  ^leed.  turn  ime.  depth,  position,  and  sea  dqitii  for  fine  tunii«  purposes,  and  finally.  (4) 
exhaustive  testing  against  as  mudi  real  data  as  possible  to  verify  or  disprove  Dwyer’s 
ooi^ectuies.  qredfically.  what  signals  are  optinud  for  kunosis  as  an  effective  detection  statistic 
in  the  active  case,  and  what  kuitosis  window  lengths  and  post  processing  are  optimal  in  tiie 
passive  case.  The  utiliQr  oi  kuitosis  for  detection  would  then  have  to  be  compared  witii 
standard  second  order  methods  in  a  ROC  curve  sense  for  tpiantification  purposes. 

It  is  worth  noting  ttiat  work  has  begun  on  an  algoihfam.  Induced  Kurtosis  Statistic 
Anomaly  Detector  (KSAD),^^  to  test  Dwyer’s  conjectures  with  die  initial  results  being  very 
encouraging.  The  algoridan,  currently  in  devdopmem  is  a  sUding  windowed  kuitosis  in  time 
over  an  acoustic  cytic.  Tests  were  performed  usiiig  data  generated  by  the  Target 

Echo,  Noise,  and  Reverberation  (TENAR)  Modd.^  The  results,  fi>r  a  wide  selection  of 
platfimn  and  taiget  dytusnic  parameters  in  shallow  water  scenarios  at  varying  taiget  aqxcts, 
indicate  successful  target  location  as  an  extrema  over  the  acoustic  cycle  in  almost  80%  of  the 
cases  widioot  any  type  of  signal  post  processing  for  a  comidete  range  of  SNRi »  IS  dB  down 
to  ‘IS  dB  in  decrements  of  5  dB.  The  algorithm  uses  no  FFTs  and,  as  sudu  is  indqiendent  of 
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dossier.  It  neidy  tties  a  nomudized  sin^ie  foimh  Older  momeit  Oauto^)  as  a  d^ect^ 
staiiaric  lo  exploit  the  active  si|^  peftttitMik»  upon  badacanetiog  as  a  function  of  target 
dynamics  and/or  modulation  effects  suggested  by  Dwyer.  In  shon,  it  is  conconed  only  with 
the  iBMisrical  distribution  of  the  echoed  signal  over  the  target  extent  Additional  shnuhtted 
tests  against  specific  taigett  of  huerest  at  aspectt  of  inteieit  for  bod)  diallow  w«er  and  deq> 
oce»  acenasios  were  petfoimed  with  a  Ugh  d^iee  of  aucoeas.  upward  of  90%.  in  the  SNRi 
range  (15.-1S]  dB.  All  simulated  tests  utilized  a  IS  degree  beam  at  a  look  angle  of  0 
degrees.  i.e..  diiecdy  on  the  Main  response  Axis  (MRA).  Preliminary  tests  against  real  data 
(static  and  dynamic)  have  also  been  peifmmed  witii  die  algorithm  again  achieving  success  in 
the  SNRi »  [IS.- 15]  dB  range  for  abbreviated  teat  sets.  More  teal  data  needs  to  be  analyzed  to 
ttioioughly  exercise  IKSAD  over  a  trial  set  oompiising  the  most  probable  scenarios  which  wU 
be  enootmiered.  TUs  is  necessary  to  complete  the  evaluation  of  Dwyer’s  conjectures. 

The  TENAR  model,^  used  for  the  generation  of  simulated  data,  predicts  the  levels 
and  complex  cross^ondations  of  sound  "heard"  by  a  roultidiannel  active  or  passive  sonar 
system  in  a  qiKcified  environment.  Sound  sources  may  indude  ambient  noise,  target  radiated 
noise,  sonar  sdf  noise,  target  echoes  and  teveibetttion  (surfooe,  vdurne,  and  bottom).  Two 
propagation  models  are  o£fored:  isovelocity  (strait  line)  wifo  boundary  reflectioos,  or 
Weirfoeig’s  CONGRATS  rqr  tiacing.^^  Eggan  and  Goddard  developed  the  sequence  of 
piQgrioa  for  generating  TENAR's  simulated  multichannd  reveiberaiion.^  These  programs 
are  fiirdier  described  in  broad  oudine  foirn  by  Luby^^  with  the  physical  and  mafltematical 
badegroond  given  by  Piinodiouse.^ 

TUs  reveibeiation  capability  has  now  been  integrared  irtto  The  Sonar  Simulation 
Toolset  (SST)^  developed  with  qponsordiip  from  sevmal  U.  S.  Navy  sources  including  the 
^ipUed  Research  Labmatory  of  The  Pennsylvania  State  UdveisiQr  (Lee  Oliver.  Leon  Sibul. 
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and  Rank  Symons.  Jr.).  SST  is  a  set  of  oxnpuier  prognnis.  input  files,  (^ject-orimiied 
software  components,  and  software  development  tools  for  building  and  tunning  sonar 
simulations.  In  this  context,  a  "sonar  simulatioo"  is  a  OMnpuier-based  process  for  predicting 
the  respcmse  of  a  sonar  system  to  a  paiticular  environmett.  It  produces  a  digital  representation 
of  the  predicted  signal  in  cadi  diannel  of  the  sonar  receiver’s  processing  padt  This  signal 
indudes  random  fluctuatkms  with  the  cmrect  statistical  properties.  The  SST  enables  a  user  to 
create  an  "aitifidal  ocean"  for  (1)  testing  new  or  proposed  sonar  systems  or  tactics,  (2) 
training  sonar  operators,  (3)  planning  eiqieiiffients,  or  (4)  validating  models  of  underwater 
acoustic  phenomena  by  comparing  simulatioo  results  with  measoiemems.  Inputs  to  a  sonar 
dmuladon  indude:  (1)  the  characteristics  of  the  ocean  itself  (sound  propagation,  sound 
scattering,  ambient  noise),  (2)  active  targets  (scattom  of  an  active  sonar’s  pulses).  (3)  passive 
targets  (noise  sources),  (4)  the  sonar  transmitter  (pulses,  beams,  trajectory),  and  (5)  the  sonar 
receiver  (beams,  trajectory,  signal  processing,  self  noise).  SST  commands  available  now  can 
generate  reverberation  and  target  edioes  for  a  multichannd  active  sonar  and  compute  the 
signal  received  by  a  multichannel  passive  srniar  from  any  number  of  maneuvering  sources. 
Three  different  methods  for  generating  reverberation  are:  (1)  tandcmi  point  scatterer  based.  (2) 
Gaussian  integration,  or  (3)  Monte  Carlo  integration.  A  suaisJit  line  Gsovelodty)  sound 
propagation  modd  is  built  in;  alternatively,  the  user  may  diorae  a  propagation  modd  based  on 
eigentay  files  produced  by  NUWC’s  (Generic  Sonar  Modd  (GSM).  The  latter  medianism 
makes  all  five  of  GSM’s  i»opagation  models  available  to  the  SST.  The  scattering  function 
produced  using  foe  Gaussian  and  Monte  Carlo  integratimis  is  a  time  dependent  cross  spectral 
dend^  matrix  for  reverberation.  SST  uses  fins  scattering  function  to  (xoduce  a  multidiannd 
stodiasdc  time  series  for  subsequent  analysis.^^ 
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It  if  the  author's  intent  to  use  SST,  in  panictilar  SST’s  ability  to  leflect  the  oonrea 
statistical  propeities  of  the  received  signal  for  a  specified  cnviroonent.  to  quantify  his  woik 
widi  KSAD^^  in  relation  to  cunent  detection  methods.  The  final  result  will  be  a  qualification 
of  Dwyer's  conjectures  regarding  the  utilify  of  time  domain  kuitosis  as  an  effective  detection 
stmistic  foe  both  active  and  passive  sonars. 
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SOFTWARE  DESCRIPTION 
A  brief  descriprioo  of  Hi-Spec  Routines  follows^ 

Hiafaer-Oider  Spectnm  Esitention:  COnventiopal  Methods 
cuniLest  Computes  sanipk  estimNes  ttf  cmiulsnts 
bispec_d  Direct  method  fi>r  biqrectnim  estimation 
biqtecj  Indirea  inethod  ibr  bispectium  estinurion 
^_stat  Detectkm  «“rt«rira  for  Hbddi’s  Oaussbrnity  and  Linearity  tests 

Hjgher-Ofder  Soectniin  Estimirin^  Pafametrio  Methoda 
cumjtnie  Confutes  theoretical  (tnie)  cumulsnts  of  AR,  MA  and  ARMA  procesaes 
fp_iid  Geneiaies  sequence  of  i  J.d.  nndom  variables 

amuLsyn  Oenmaes  ARMA  synthetics 

ar_rcest  AR  paiatneteis  ooirelation  andAor  «wniii*my 

tnajcst  Estimates  MA  paiatneteis 

ann^qs  Estimates  AR  panmeien  using  the  q-slice  algoridan 

amuLits  Estimates  ARMA  paiameten  usfaig  the  residual  time  series  method 

tls  Total  least  squares  stdutioo  to  a  set  of  linear  equations 

Quadratic  Phase  Courting 

qpc_gen  Genentes  quadtmically  phase-coupled  harmonics  in  noise 

qpc_lor  Paiametiic  QPC  detection  of  quadratic  phase  coiqding  via  die  TOR  mediod 

Hannonic  Retrieval 

hami^en  Generates  hannonics  in  Gaussian  (cokxtd)  noise  for  the  harmonic  retrieval  problem 


tde_fen 

ttk 


dcMLjea 

dot 


rivjD 

riv_tr 

iv.cal 

bioq» 
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Ettimtrion  of  power  ^seoit  of  hannooks  using  MUSIC,  Eigenvector, 
Pisaicnko,  ML  (Capon)  and  AR  methods  baaed  on  the  diagonal  slice  of  the 
fooidHxder  cumulant  wldi  conventional  periodogram  also  given  ftn* 
oon^Mtison  poipoaes 


Time-Delav  E«rim«tton 
Oeneiaaes  Synthetics  for  tuncHlelay  estimation 

Estimatrs  time-delays  from  two  sensor  measmements  using  dtiid  order  cross 

ciwMilawia 


Array  Processing 
Oeneraies  Synthetics  fin  DOA  pioblem 

Estimates  number  of  sources  and  their  bearings  fiom  a  linear  array  of  sensois 
using  MUSIC,  eigenvector.  Pisarenko,  ML  (Capon)  and  AR  mediods  based  on 
qiatial  fbuitb-order  cuinulants 

^  Prediction 

Adaptive  LP  using  double  lattice  filter 
Adtqttive  LP  using  transversal  filter 
Conqaites  instiumetttal  variable  processes 


Magnitude  and  Phase  Retrieval 

Estimates  impulse  reqxnse  via  tiie  bioqntrum  method  and  cooqnites  conqilex 
cqtstrum  of  a  signal 
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ouAil  Fottfkr  Ftauc  wd  iMgMtudc  of  1  Signal  usiQg  the  Matsuoka-Ulfych 

algofitfam 


cuni2_est 

cmnSjcst 

cuiD4jest 


Undocumented  R 
Estimatef  oovarianoes 


Estimates  ddid-oider  cumulaite 
Estimates  fouith-oider  cumulnts 


A.1  BISPECTRUM  ROUTINES 

The  biqtectnim  of  the  process,  y.  is  estimated  via  the  direct  (FFT-baaed)  and  indirea 
metfaods.^*^^ 

For  the  direct  medud.  te  time-series,  y.  is  lamented  into  possibly  overlapping 
records;  tte  mean  removed  fiom  eadi  reotxd.  and  the  FFT  computed;  the  bispectnim  of  die  k- 
di  record  is  oooqMted  as 


\ 


(A.M) 


wbea  denotes  the  FFT  of  the  Adi  record.  The  bispectnl  esthnairs  are  avenged  across 

records,  and  an  optional  frecpiency-domain  smoother  (Rao-Gabr  window)^^  is  applied. 

For  the  indirect  method,  dw  time-series,  y.  is  again  segmented  into  possibly 
ovetligiplng  records;  unbiased  samide  estimates  of  third-order  cumulants  are  cougiuted  for  eadi 
record  md  dien  averaged  across  records;  a  Parzen  lag  window  is  applied  to  die  estimaed 
gnmHiwua,  and  die  bispectium  is  obtained  as  the  2-D  FFT  of  die  vdndowed  cumulant  function. 
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AJt  GAUSSIANITY  TESTS 

A  decision  statistic  for  Hinich’s  Gaussianity  test  is  estimated.^  The  bi^wctnim  is 
estimated  using  the  direct  method,  and  a  frequency-domain  2-D  boxcar  smoother  is  apfiiied. 
The  power  spectnmi  is  estimated  via  the  direct  method,  and  a  boxcar  smoother  is  ^tplied.  The 
bicoherence  is  then  estimated.  The  Gaussianity  test  (acuuUy  xeiD-skewness  test)  basically 
invtrfves  deciding  whether  or  not  the  estimated  bicoberence  is  zero. 

The  statistic  for  the  Gaustianity  test,  sg,  is  a  three-element  vector  where  sg(\)  is  the 
estimated  statistic  S.  sgO),  is  the  number  of  degrees  of  freedom,  p.  and  the  probability  of  fitise 

tdaim.  SgO),  is  the  probability  that  a  ^  random  variable  with  p  degrees  of  freedom  could 

have  a  vidue  laiger  than  the  estimated  5  in  sg(lX  If  tttis  probability  is  smdl,  say.  0.05.  then 
we  may  reject  the  hypothesis  ot  zero  skewness  at  a  FFA  (or  level)  of  0.05. 

Therefore,  if  you  decide  to  accqpt  the  hypothesis  that  die  data  have  non-zero  skewness,  then 
the  probability  diat  the  data  may  actually  have  zero  skewness  is  givca  by  sgO\  If  FFA  is 
large,  then  dre  hypodiesis  of  zero  skewness  cannot  be  easily  ityecred. 

The  normalized  bispectram  (bicoherence)  is  defined  as 

- —  (A2-1) 

(S(«,)S(»^«, ♦«,))«* 

^  ^  bispectrum  and  S(q)  is  the  power  spcctnim,  Under  the  Gaussianity 

(zero  skewness)  assumption,  the  e*^  'ted  value  of  Is  zero.  The  test  of  Gaussianity 

is  based  on  die  mean  biooherenoe  power. 


where  5(fc)  is  the  Kionecker  deha  fiinctiofL 
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The  jUt  power  spectrum  is  obtained  as  follows:  First,  a  rank  p  ^)proximation  of  the  matrix  C 
is  obudned.  as  ^  ■  ySv\  where  ^  is  obtained  from  S  by  setting  «0,  k  « 

The  AR  parameter  vector  is  then  obtained  as  the  stdution  to  (2a  >  Q;  die  method  in 
Cadzow^*  is  used,  and  the  solution  is  forced  to  have  unity  modulus.* 

The  ML(Capon)  solution  is  given  by. 


(A.3-3) 


6] 


INTRODUCTION 

This  qjpendix  presents  an  overview  of  the  definitions,  prqpeities.  and  estimation 
methods  of  higher-oider  statistics  or  cumuhmts  of  stochastic  signals  nd  their  associated 
Fbuiier  tiaiafoims  known  as  higher-order  q)ectni  or  poly^KCtra  as  presented  in  the  tutorial 
refisrences.^*^*^® 


B  J  MOMENTS  AND  CUMULANTS 

Higfier-otder  spectra  of  stationary  stochastic  signals  are  defined  in  terms  of  cumulants 
and  are  referred  to  as  cumulant  spectra.  Given  a  set  of  a  real  random  variables 

their  joint  cumulants  of  order  r  •  *  ...  *  defined  as 


(B.1-1) 


S 


0, 


vidiere 


•  (OpUj,  o J  »  E  iexpjio^  +  ...  (B.1-2) 

is  their  joint  characteristic  function.  The  joint  moments  of  order  r  of  the  same  set  of  random 
vaiiaUes  are  given  by 


i-ff 


d^^(U|y  ^2* 


dJ} 


(B.1-3) 
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Tbe  joiitt  cumulants  can  be  exfMcssed  in  tenns  of  the  joint  moments  of  the  set  of  random 
variables  as  follows  for  the  moments 

=  Etcj)  iNj  =  £{xi*} 
nij  «  £{xf)  «  £{x^ 

of  die  random  variatde  {x^}  are  related  to  its  cumulants  by 

c^  •«, 

Cj  -  -  mf 

e,  «  m,  -  SnijiR,  *  Im* 

C4  •  -  3m^  +  I2m^l  -  Sm* . 

It  is  clear  diat  die  computation  of  the  rth-order  cumulant  cf  (x^)  requires  knowledge  of  all  its 

moments  from  order  first  to  rdi,  i.e.,  k  *  0,  ±1.  ±2, ...  is  a  real, 

strictly  stationary  random  process  with  zero  mean.  E{x(k))  s  0,  then  the  moment  sequences  of 
the  process  are  related  to  its  cumulants  as  follows: 


(B.1-4) 
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wdieie  the  moments  up  to  onkr  n  exist  and  depend  only  on  the  time  differences 

Tj,  T2,  _  |.  Taking  into  account  the  zero  mean  condition  ^x>ve.  the  cumuhnts  are 

related  to  the  moments  as  follows  for  mders  n  •  123.4: 

Cl  -  Ml  -  £{x(*)}  -  0  (msan  vabu),  (B 

c/(t,)  »  m/(T,)  -  =  »i2*(T,)  »  E{xik)xik  +  t,)}  ^  ^ 

(covariance  sequMiee), 

udiere  (f  is  die  autocoireladon  sequence. 

«s*(^i.^a)  •  "  'iM  *  2m,** 

(B.1-5C) 

(Aird--order  rnomemoi  cumukaa  sequence), 

<^4  (“ft. ‘'2*^3)  =  »^4*(V‘f2»^3>  -"b'C-Cl)  •"2*(^s  -  *^2)  ■  "»2  (“^2) 

•m2*(T,  -  T,)  -mj'Ctj)  -m/Ctj  -  t,)  -  m,*[m3*(T2  -  Tj.Tj  -  t,) 

♦  m/Ct^.T,)  ♦  m/(x2,t^  ♦  «3*(^i»^2)J 

♦  2m,**  [m/(T,)  +  m/Cx,)  ♦  m/Cx,)  +  nij'Cx,  -  x,) 

+  m/ (x,  -  +  m^ix^  +  t,)3  -  6m,*^- 


(B.l-5d) 
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^2»  ^3)  "  ***4  ^9)  '*  **2  (^1)  *  ^  (^3  ”■  Tj)  -  III2  (Tj) 

•  "  T,)  -  Wj'Cts)  ■  "fi) 

»  E{xik)x(,k  *  T,)x(ik  +  Tj)} 

(fixirth-order  moment  sequence). 

Equations  (B.l-5aHB.l-5d)  illustrate  tfiat  while  the  mnnents  up  to  third  order  are  identical  to 
cumulants  for  zero  mean  stationary  random  processes,  the  generation  of  the  fouittw>rder 
cumulant  sequence  requires  knowledge  of  the  fourth-order  moment  and  autocorrelation 
serpienoes. 

By  putting  the  delays  3  0  in  equations  (B.l-5aHB.l-5d)  with  the  zero  mean 
assumpdon,  the  variance,  skewness,  and  loutosis  becmne 

Yz  *  <2  (0)  *  (yoriance) 

y5  •  ciaXO  -  £{x*(ik))  (sfownesu)  (B.1-6) 

yJ  -  C4*(ao.Q)  -£U^(*))  -  3[y5*  ikurtosis) 

B2  HIGHER-ORDER  SPECTRA 

Ifigher-order  spectra  of  stochastic  signals  are  usually  defined  in  terms  of  cumulants 
and  not  mranents  for  two  reasons:  (1)  for  a  statitmary  Gaussian  random  process,  all  of  the 
moments  for  n  >  2,  although  generally  non-zero,  provide  no  new  informaticm,  while  the  fact 
that  the  joint  cumifiatus  for  n  >  2  are  idenfically  zero  conveys  tlus  explicitly;  (2)  if  the  random 

variables  {X|,...,x^}  can  be  divided  into  any  two  or  more  groups  udiich  are  statistically 


6S 

indcpeadent,  tbdr  nth-onkr  cumulants  tit  identictUy  zero  ttaus  providiiig  t  metsuie  of 
stadstictl  independence.  (The  nth-onkr  moments  ait  non-zero.) 

Assuming  that  the  cumulant  sequence  is  bounded.  i.e.,  satisfies  the  condition 

E  •••  E  . V,)|<%*>*n»»w>«tercumulantspectnimc;(o,,...,«^.^^ 

—  Vi  ■  - 

of  (xfA)}  exists  and  is  continuous  and  is  defined  as  the  (n-l)-<limensional  Fourier  transfoim  of 
the  nth-onkr  cumulant  sequence 


C’li  (®i»  •••» 


1 

(2irr‘ 


•exp (-/(«! tj  *  ...  + 

I  i  >1,2,...,r-1  and  |0|  +  O2 


In  general.  C/(0|,  ...» o..i)  is  complex  for  n  >  2.  ix..  it  has  magnitude  and  phase. 


®«-i)  "  I  C®!***’**"**-!^  I  til  (®|f  •••»  *^*-1)  • 


(B.2-2) 


The  cumulant  spectnim  is  also  periodic  with  period  2n.  i.e.. 


C/Co, 4  2n,...,«,_|  ♦  2ii). 


The  power  spectnim.  Inqiectnim.  and  trispectxum  ait  special  cases  of  die  nth-onkr  cumulant 
qwctium  defined  by  equation  (B.2-1). 
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Pvwtr  SjMdnuH.*  mm2 

C/  («)  «  ^  Cj*  (t)  exp  { -J  («t)  }  (B.2-3) 

|o|  ^1(. 

where  c/(t)  is  the  oovarince  sequence  oi  {x(k)]  given  by  equ^on  (B.l-5b).  From  equation 
(B.l-Sb)  and  equation  (B.2-3)  we  have 

Ca*(T)  -  c/(-T) 

C/(«)  -  Ci*(-«) 

>  0  (leil,  moimegative  fiaetiOH). 

Bitpmamm:  mmj 

Cs* («i. «a)  =  — ^ 

«i  • 

|«,  +  Ojlix, 

(B.2.5) 

where  t,)  >>  third-order  cumulam  wquence  of  lx(k)}  described  by  equation  (B.l- 


5c).  Iiiqx)rtam  symmetiy  conditions  follow  fmn  equation  (B.l-Sc) 
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(^i»  ^  ^  "  ^a) 

■  Cs*(t,  -  T,,  -  T,)  -  C5*(t,  -  tj*  “  ^a) 

«  -  T,). 

Tlw  definition  of  the  bupectnim  in  equation  (B^-5)  nd  the  ptopeities  of  thitd-oider 
cunulants  in  equation  (B  J-6)  give 

C/(tt|.Qa)  -  C,'(02.0|)  -  (-«y  -«,) 

»  “»j)  ■  C5*(-«,  -  »r^a) 

(B.2-7) 

«  -  «a»®i> 

-  C/(Oy  -  «j). 

TUa  knoudedge  of  the  biqKcnum  in  the  triangular  region  k  0,  Uj  k  Uj*  ^a  ^  ^ 

enough  for  i  contpleiB  description  of  the  bispecminL 
Trispectnun:  h»  4 


49  49  49 

c/(», ^  E  E  E 

C2»)*  T|-  —  «,•  —  »*•-• 

+  «,  +  «,  1^*, 


(B^-8) 


where  c/  (T|,  t,,  )  >*  foe  fouitlKHder  cumulam  sequence  given  by  equation  (B.l-Sd). 

Symmeuy  properties  can  be  derived  for  foe  triqtectnun,  rimilar  to  those  given  in  equation 
(BJ2-7)  for  the  biqrectiuin. 
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A  normalized  cumulam  spectnim  which  ccxnbines  the  cumuUnt  qiectium  of  order 
II,  C/(U|, power  spectnim,  C,*(<i>).  of  a  process  may  be  defined  as 


wt*  *  \  ^*^1* 

(«i.  ; - - - ;; - ; - 


(B^-9) 


The  thiid-order  (n  s  3)  case  is  called  die  bicoheience  which  is  the  normalized 


Uqiectnim 


_ Q*  (0|,  tta) _ 


(B.2-10) 


wdiere  C/  is  defined  in  equation  (B,2-S)  nd  is  defined  in  equation  (B2>3). 


The  six  symmeby  regions  of  the  third-order  cumulant  as  given  in  equation  (B2-6)  and 
die  twelve  symmetry  regions  of  the  biqiectrum  as  given  in  equation  (B2-7)  are  listed  below 
and  shown  in  Figures  B  J-la  and  B2-lb: 
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cumulnt  symaMries 

5.  Cj*(T,-T,.-T,) 
6e,*(T,-Ty-Tj) 


bupecmm  synmieines 

*•  ^*(U|,  O3) 

3.  Cj*  (-«y  0|  ♦  «j) 

Cj*  (-«,, «,  ♦  «j) 

5.  Cj*(-«,  -U2>»t) 

6.  C^* (-«,  -  Oy 

7. C^(«,.-03) 

8-  (-Wp  -«,) 

9.  C,* (wj,  -  «j) 

10.  c/  (Oj,  -o,  “  «j) 

11.  +Oj,-Oj) 

12.  (0|  ♦  U2»  -Wj) 
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Figure  B  J*le  Figure  B^lb 

The  dx  qruunetrj  regions  of  the  The  twehre  qrnunetry  regions  of  the 

third  order  cumulant  ftmction  e|*(t|,  tj)  Uspectrun 

The  additional  synunetiies  as  defined  by  the  extiemhies  of  the  outer  hexagon  in  the  above 
figure  are  due  10  die  effects  of  a  discrete  time/continuous  fretpiency  biqiectnim  as  opposed  to 
die  standard  inner  hexagon  symmetiy  r^on  of  the  continuous  time  bi^Kctnun  oorreqxmding 
to  a  band  limited  process.  The  net  result  is  an  additional  set  of  twelve  symmetry  regions 
defining  the  borders  of  the  outer  hexagon.  The  douUe  periodicity  of  the  biq)ectrum  causes 
this  effect  and  is  illustrated  in  a  Figure  B.2-2  fiom  Allisoa^  The  Uodt  cr^tal  "C  notation 
defines  conjugate  symmetry  regitms  of  the  biq^ectrum  vdiOe  the  boxed  numbers  donate 
etpnvalem  regions  of  the  discrete  time  Ix^wctrum.  Hgure  B2-2  reflects  a  normalized 
samitiing  fiequency  •  1  with  the  lighter  shaded  region  being  the  zero  part  of  the  bispectrum 
for  a  bandlimited  luooess  indicative  of  a  properly  sampled,  stationary  time  series.  The  darker 
diaded  region  refuesents  die  non-zero  part  of  the  InqiecttuiiL 


Figure  BJ*2  Fundamentil  Domiu  of  Dbcrelc  Time  Biipeclnun 


Previous  U^pectnim  properties  ^){dy  solely  to  ieal>vdued  signals.  The  extension  to 
comidex-valued  signals  has  not  been  widely  addressed.  The  following  definitions  and 
properties  are  due  to  the  requirement  for  cmnidex  signals  for  radar  signal  processing  as 
presented  in  Jouny.^  This  is  a  normal  lequiiemem  finr  analytic  signal  analysis  for  a  wide 
qpectnim  of  apfdicatimis  involving  real  or  quadreture-sami^  data. 

If  the  data  {x(ik)}  is  a  comidex  stationary  process,  the  seemtd  order  cumulam  is  defined 

**  ”  E{x*(k)x0c  *  t)}.  is  also  possible  to  ahemate  the  conjugate  so  that 

c/ (t)  ■  E{xik)x*(,k  *  t)  }.  third  joim  cumulant  case,  the  conjugate  can  be  placed 


either  on  one  or  two  entries  of  die  triide  inodua  in  equadrm  (B.l-5c).  Only  one  placement. 
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however,  defines  which  one  of  the  symmetry  leUdoos  given  in  equation  (B.2'7)  remains  valid. 
In  general,  the  bispectrum  of  complex  signals  has  twofold  symmetry  about  an  axis  defined  on 

»  (02><i>2  “  •  -.SOyO,  ■  -Oj).  See  Rgure  B.2-3  for  the  distributions 

of  responses  in  the  bispectral  domain  as  a  function  of  the  (Aacement  of  the  conjugates  for  the 
third  order  cumulants. 


«2 


Figure  B  Distribution  of  the  responses  in  the  hispectrai  domain  as  a  fiinction  of 
third  order  cumuiant  cortiugate  piaoement 

For  example,  if  the  first  quadrant  ^  preferred  then  the  third  order  cumuiant  is 

defined  as  Tj)  *  Eix*  ik)x(k  +  x^x(k  *  tj)). 

The  following  relations  define  the  bispectrum  of  complex-valued  signals  for  all 
possible  positions  of  the  complex  conjugate: 


With  the  symmetry  reUtioa  C/(«,,Oj)  -  C/(«2, «,), 

£(x(l}x'(**t,)x(*»t^)  ^.5.,  <*’(-o^X(o^X(-u,-o^> 

With  the  symmetry  lelatkm  -Wj  -  Wj), 

£{i«Qx«  ♦  t,)x*  (Ik  ♦  x,)l<  <  I(opX'(-u^X(-o,  -  > 

With  the  symmetry  idatian  C/COj.Oj)  - 

Eix^x^Qt  * x^)x•(fe  +  }< .f.> <X-(-o,)X’(-«3)X(-«,  -  «,) > 

with  the  symmetry  rdatioa  -  C/(02»UiX 

£lx*(*)x*(lk  ♦  T,)x(t  ♦  .f.»  <f‘  (-«P  X(«j)  X*  («,  »  > 

with  the  symmetry  relation  C^(o,,«^  «  €,*(-«,  - 

£(x’(t)x(k  ♦  T,)x*(t  ♦  tp) ,  .f_,<X(o,)X’(-»,>X*(«,  t  u^> 
with  the  symmetry  relation  Cj'COpOj)  » 

F 

< - > 


«4iete 


denotes  die  Fourier  transftxm  pair  and  o  denotes  die  ensemble  average. 
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The  triple  ooiieUtion  aq)ect  of  the  iMspectiuin  is  self  evident  fiom  its  definition  as  the 
Fburier  tiansfonn  of  a  third  oner  cumulant  in  eqiuoion  (B.2-1).  This  definition  ini{dies  that  the 

bispectnun  at  a  patticular  frequency  pair  (u^up  is  non-zero  if  at  least  one  of  the  three  q)ectral 
I'c^’onaes  ^  up  >4  correlated  with  the  other  two  since 

Cj^Cca^up  ■  <jr(QpX(upX*(Uj  +  up>.  This  definiticm  arises  ftom  a  Fburier-Stieltjes 

representation  of  X  where  the  distinction  is  made  between  the  power  ^tectnun  iqmsentation  of  the 
contributkn  to  the  mean  product  of  two  Fourier  components  whose  frequencies  are  the  same 
conqMied  to  the  similar  bispectnun  contribution  invr^ving  three  Fourier  components  where  one 
frequency  equals  the  sum  of  die  other  two.^^ 
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APPENDIX  C. 

SIMULATED/REAL  DATA  TEST  RESULTS 
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POWER  SPECTRUM/BISPECTRUM  COMPARISON  FOR  SIMULATED  DATA 
Figures  C.0-1  thiDU^  C.0-3  show  a  comparison  between  the  power  spectrum  and 
bispectrum  of  a  harmonic  (pure  tone)  for  three  different  states:  (1)  noiseless  (Figure  C.D-l),  (2) 
low  noise  (Figure  C.0-2),  and  (3)  high  noise  (Figure  C.0-3).  The  low  and  high  noise  stales 
represent  signal-to-noise  ratios  (SNR)  of  16.S  and  *3.5  dB  respectively  where  the  SNR  is 
defined  as  the  ratio  of  the  signal  variance  to  the  noise  variance;  when  expressed  in  dB,  it  is 

given  by  10  log  (o^o^)  where  af  is  the  variance  of  the  signal  and  is  the  variance  of 

the  noise.  The  bispectrurri  was  computed  using  a  FFT  size  of  256  corresponding  to  the  data 
length  as  no  averaging  was  performed.  The  bispectrum  slices  reflect  the  effects  of  no 
windowing  aixl  a  Rao  Gatv  window  of  length  S  for  smoothmg  purposes.  The  slices  with  the 
highest  bispectral  value  in  the  triangular  region  of  support  were  chosen,  a  convention  followed 
in  bi^iectral  procesang.^  The  pure  tone  was  designed  to  show  up  at  .5  hz  for  this  set  of  test 
cases  with  a  samiding  rate  of  2  hz.  The  noise  states  were  generated  by  adding  Gaussian  noise 
to  the  pure  tone.  The  plots  indicate  essentially  no  difference  in  the  signal  peak  to  noise  level 
representations  between  the  power  spectrum  and  bispectrum  for  the  noise  states  in  Figures  C.0- 
2  and  C.0-3  reflectively. 

C.1  DIRECT/INDIRECT  BISPECTRUM  COMPARISON  FOR  SIMULATED  DATA 
Hgures  C.l~l  and  C.1-2  show  a  comparison  between  the  conventional  direct  and 
indirect  mediods  for  estimating  the  bispectrum  as  defined  in  sections  2.1.1  and  2.12 
respectively.  The  data  for  the  figures  was  generated  using  the  Hi-Spec  program  "qpc.gen” 
i^ch  generates  quadradcally  phase-coupled  harmonics  in  noise  as  defined  in  section  A.O  of 
Appendix  A.  Hgure  C.1-1  rfnesents  laspectial  processing  using  a  quadradcally  phase  coupled 
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synttietic  oonsistiiig  of  ttiree  quadnticaUy  phase  oouf^ed  haimonics  at  fiequencies  fl  »  0.12.  f2 
»  0.18,  and  f3  s  n  t2  «  0.30  hz.  The  data  was  segmetued  to  reflect  independent 
realizations.  FOr  eadi  realizatkm  the  phases  of  the  first  two  haimonics  were  chosen  randomly 
with  the  phase  of  the  third  haimonic  set  to  the  sum  of  the  phases  of  the  first  two.  A  total  of 
64  independent  realizations,  each  consisting  of  64  samples,  were  generated.  The  amplitudes  of 
aU  three  harmonics  were  unity  and  the  signal  was  noise  free.  Bgure  C.1-2  data  was  idemical 
to  Figure  C.1-1  with  the  exception  that  Gausdan  noise  was  added  to  the  signal  to  obtain  an 
overall  signal-to*noise  ratio  of  0  dB. 

The  direct  bispectrums  from  left  to  rigltt  reflect  no  windowing  and  a  Rao^labr 
window  of  length  5  respectively.  Similarly,  the  indirea  iMSpectiums  from  left  to  right  reflect 
i»  windowing  and  a  Paizen  window  respectively.  The  FFT  size  in  both  cases  was  128  with 
the  number  of  lags  chosen  in  the  indirea  case  for  die  estimation  of  the  higher  order  moments 
being  21.  (It  is  suggested  as  a  rule  of  thumb  in  Hi-^iec  that  die  number  of  lags  should  be  sa 
to  the  number  of  samfdes  per  segment  divided  by  10.) 

The  severe  deterioration  of  the  indirea  Uqiectnim  estimate  for  the  noisy  data  in 
Rgure  C.1-2  was  unexpected  reladve  to  the  successful  suppression  of  noise  noted  for  the 
direa  bispectium  estimate  in  comparison.  The  skewness  of  the  noiseless  data  in  Figure  €.1-1, 
a  funcdon  of  the  quadratic  phase  oouiding,  was  0.82  indicating  highly  skewed  data.  The 
skewness  of  the  noisy  data  in  Figure  C.1-2  was  0.36  indicadve  diat  the  phase  coupling 
contribution  to  die  skewness  was  na  masked  by  the  addition  of  Gaussian  noise.  It  is 
worthwhile  to  note  diat  for  the  generation  of  uncoufded  situisoids  at  the  same  frequencies. 
0.12, 0.18,  and  0.30  hz,  the  noiseless  skewness  was  0.03  while  the  noisy  skewness  was  -0.01 
indicative  of  symmetry  whidi  exfdains  die  inability  of  the  bispectrum  to  differentiate 
unooiqded  sinusoids  from  symmetric  pdf  noise. 


78 

The  indirect  bispectnim’s  inability  to  sufficiently  suppress  the  noise  in  Figure  C.1-2 
along  with  resdution  proUems  noted  in  Figure  C.1-1  (note  the  loss  of  the  quadratic  phase 
coupled  triplet  in  the  Parzen  window  case)  dictated  the  use  of  the  direct  bispectnun  for 
comparison  against  second  order  statistical  methods  for  this  pqier.  The  proUems  related  to 
the  indirect  bispectium  remain  unresdved  at  this  time  pending  further  study  in  regard  to  the 
theory  and  algorithmic  imfdementation  in  Hi-Spec. 

CJ,  POWER  SPECTRUM/BISPECTRUM  COMPARISON  FOR  REAL  DATA 

Figures  C^-1  and  C2-2  show  a  cmnparison  between  the  power  spectrum  and  direa 
biqwctrum  as  defined  in  equations  (B2‘3)  through  (B2-7)  respectively  of  Appendix  B.  The 
biqrectrum  was  computed  using  a  FFT  size  «  128  with  a  Rao-Gabr  window  of  length  3. 

Figure  C^-1  is  a  sonar  ping  using  a  1(X)  ms  hop  code  designated  Ping  A  while  Figure  C.2-2  is 
a  11  ms  pure  tone  designated  Ping  B.  The  1-quadrant  bi^rectrum  is  shown  via  contour  {dots 
with  the  quadrant  symmetry  clearly  evideitt  as  illusttated  by  the  figures  in  section  B.2,  in 
particular  Figure  B2-3.  "Distribution  of  the  responses  in  the  bispectral  domain  as  a  function  of 
third  order  cumulant  conjugate  idacement."  udiidi  is  in  the  first  quadrant  for  Hi-Spec. 
Bispectrum  slices  are  computed  following  die  convention  referenced  in  section  C.0  for 
comparison  with  the  power  ^rectrum.  The  ping  data  limitations  constrain  the  averaging  to 
approximately  9  segments  minors  die  tests  performed  in  Wilson.^  The  frequency 
scales  have  been  normalized  for  sanitizatitm  purposes  and  the  power  qiectrutn/bispectrum 
comparison  reflects  a  relative  dB  scale  as  a  log  power  spectrum/log  magnitude  bispectnim  is 
computed  -  a  parameter  typically  found  in  the  literature. 

The  figures  indicate  that  no  imxxssing  gain  is  realized  by  the  biqiectrum  over  the 
power  spectrum  for  detection  purposes  in  both  instances.  Goser  examiruttion  of  the  statistics 
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fcveals  dutt  die  distribudons  of  the  letums  had  dcewness  measures  of  (-0.02. 0.00)  and  (0.02, 
0.01)  respectively  for  their  real/imaginary  parts.  Moreover,  the  kuitosis  statistics  were  (3.68. 
3.1S)  and  (3.19.  3.21)  respectively  (a  kuitosis  of  3  indicates  Gaussianity).  Oeaily.  the  desired 
effects  of  taiget  dynamics  and/or  modulation  necessary  to  skew  the  dam  suffidentiy  for 
bispectral  enlunoement  via  suppression  of  Gaussian  or  non<iaussian  symmetric  pdf  noise  were 
not  enough^"^^®. 

The  dam  analyzed  represented  shallow  water  scenarios.  A  total  of  approximatdy  1 1 
pings  were  with  the  results  in  ad  cases  similar  to  dmse  described  in  Hgures  C.2-1 

and  C.2-2.  The  average  echo  skewness  measure  over  1 1  pings  was  (0.01.  -0.00)  widi  a 
ricewneas  range  of  (1-0.03  0.03].  [-0.03  0.01])  and  a  standard  deviation  of  (0.02. 0.01).  The 
average  echo  kuitosis  measure  was  (3.47, 3.47)  widi  a  kuitosis  range  of  ([2.67  4.84]. 

[160  4.3S])  and  a  standard  deviation  of  ^.73. 0J4).  These  numben  appear  to  be  well  witiiin 
ttie  bounds  of  Gaussianity  for  skewness  with  only  a  sU^  variance  for  kuitosis.  The 
conesponding  mimbeia  for  the  leveibeiation  data  in  the  nd^iboihood  of  the  echoes  was  0).OO. 
0.(X))  for  skewness  widi  range  ([-0.01  0.01],  [-0.00  0.01])  and  standard  deviation  (0.01, 0.00) 
which  is  die  crux  of  the  problem,  i.e.,  the  echoes  cannot  be  discriminated  suffidentiy  in  a 
statistical  sense  from  the  leveibetation  for  successful  laspectial  piocesdng. 


Power  Spectrum 
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Figure  C.0*l  Power  Speclrum/Bispecirum  Comparison  for  Simulated  Data  -  Noiseless 


Power  Spectrum 


Figure  C.0'3  Power  Specirum/Bispeclnim  Cumparisun  fur  Sinwiaied  Data  -  High  Nt>ise 


Pinq  A  Envelope  Direct  Blspectrum 


Figure  C.2-1  Power  Speclrum/Dispcclmm  Comparison  for  Real  Data  -  Hop  Code 


Pina  B  Envelope  _ _ Direct  Bispectrum 


Sprclrum/Uispectrum  Comparison  for  Real  Data  -  Pure  Tone 


